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Abstract  of  Dissertation  Presented  to  the  Graduate  Council 

in  Partial  Fulfillment  of  the  Requirements  for  the  Degree  of 

Doctor  of  Philosophy 

SPATIALLY  DEPEIJDENT  TRANSFER  FUNCTION 
OF  NUCLEAR  SYSTEMS  USING  CROSS  CORREUTION  LETHODS 

By 

Chester  D.  Kylstra 

December  21,  1963 

Chairman:  Dr.  Robert  E.  Uhrig 

Major  Department:  Nuclear  Engineering 

The  concept  of  a  transfer  function  for  a  nuclear  sys- 
tem is  extended  to  include  spatial  effects.  The  general 
equation  is  derived  using  the  time -dependent  Fermi  age 
and  diffusion  theories  for  a  single-region,  isotropic, 
homogeneous  medium.  The  fluctuation  of  the  thermal  neu- 
tron density  at  any  point  in  the  assembly  is  related  to 
the  variation  of  the  fast  neutron  source  intensity. 

The  general  transfer  function  equation  is  specialized 
for  several  cases,  including  the  case  of  a  point  source  in 
a  cylindrical  medium.  Theoretical  curves  are  calculated 
for  multiplying  and  non -multi plying  media  and  compared  with 
the  commonly  used  lumped  parameter  transfer  function. 

The  results  indicate,  in  general,  that  the  lumped 

parameter  model  predicts  the  correct  behavior  of  the 

nuclear  system  only  if  the  output  detector  is  located  a 

specific  distance  from  the  source.  If  the  detector  is 

located  elsewhere,  the  lumped  parameter  model  is  not 

capable  of  accurate  results. 
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Experiments  were  performed  on  light  and  heavy  water 
subcritical  assemblies  to  measure  the  spatially  dependent 
transfer  function  between  two  detectors.  The  input  to  the 
assemblies  was  provided  by  a  neutron  generator  which  was 
turned  on  and  off  in  a  pseudo-random  manner  by  the  maxi- 
mum length  output  sequence  generated  by  a  shift  register. 
Two  neutron  detection  systems  provided  analog  output 
voltages  proportional  to  the  neutron  density  fluctuations 
at  different  locations  in  the  assembly.  The  voltages  were 
sampled  at  a  high  rate  by  a  data  acquisition  system,  which 
stored  the  data  on  magnetic  tape.  The  data  was  analyzed 
by  cross  correlating  the  two  outputs  and  calculating  the 
cross  power  spectra. 

Statistical  equations  are  developed  which  show  that 
the  transfer  function  may  be  obtained  by  computing  the 
ratio  of  the  cross  power  spectra  for  two  different  runs, 
with  one  detector  moved  to  a  new  location  for  the  second 
run.   This  method  of  statistical  analysis  has  several 
advantages.  Among  them  is  the  ability  to  eliminate  the 
effects  of  the  source  spectra  and  the  effects  of  the  mea- 
suring system  from  the  results. 

A  comparison  of  the  experimental  and  theoretical 
transfer  functions  provided  inconclusive  evidence  as  to 
the  validity  of  the  theoretical  model.  Other  effects,  such 
as  reflection  of  the  neutrons  from  the  surroundings,  pre- 
vented complete  agreement.  The  equations  used  to  predict 
the  behavior  of  a  light  water  system  clearly  must  take 


into  account  the  effects  of  a  distributed  source,  rather 
than  a  point  source.  The  heavy  water  data  was  affected 
by  the  small  size  of  the  system  and  the  nearness  of  other 
structures,  but  the  theoretical  curves  still  agreed  very 
well  with  the  experimental  data. 
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CHAPTER  I 
INTRODUCTION 

The  principal  purpose  of  this  study  is  to  develop  a 
Spatially  dependent  transfer  function  for  nuclear  systems, 
using  the  time  dependent  Fermi  age  and  diffusion  theories. 
The  effects  of  neutron  energy,  time,  space,  and  delayed 
neutrons  are  considered.  The  experimental  spatially  de- 
pendent transfer  function  of  a  light  water  and  a  heavy 
T.;ater  subcritical  assembly  was  measured  and  compared  with 
the  theoretical  transfer  function.  The  multiplication 
factor  of  both  systems  was  varied  from  zero  to  near  one. 

The  data  collected  for  the  measurement  of  the  ex- 
perimental transfer  function  was  processed  by  cross  cor- 
relation methods.  A  secondary  purpose  of  this  thesis  is 
to  demonstrate  the  advantages  and  accuracy  of  the  parti- 
cular method  used. 

The  transfer  function  equations  commonly  used  for  nu- 
clear reactors  are  derived  from  the  time  dependent  dif- 
fusion equation,  after  the  spatial  dependence  has  been  re- 
moved by  assuming  that  the  flux  shape  is  the  fundamental 
spatial  mode-   The  reactor  is  treated  as  a  "black  box", 
or  as  a  lumped  parameter  model  {1,    2).     This  transfer 
function  works  surprisingly  well  (J.),  especially  for  low 
frequencies  where  the  neutron  flux  appears  to  respond 
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simultaneously  throughout  the  entire  reactor  and  for  cri- 
tical systems,  where  the  band  width  is  narrow  and  the 
higher  frequencies  are  greatly  attenuated. 

A  reactor,  however,  is  not  actually  a  lumped  para- 
meter system.  The  distance  between  the  input  and  output 
devices,  as  well  as  their  relative  locations  in  the  sys- 
tem, can  have  a  large  effect  on  the  measured  attenuation 
and  phase  shift  of  a  disturbance  as  it  propagates  through 
the  system.  As  higher  frequencies  are  present  and  of  in- 
terest, as  in  neutron  wave  and  random  noise  experiments, 
the  time  and  spatially  dependent  solution  of  the  neutron 
lux  behavior  is  needed  to  adequately  describe  experimen- 
tal results.  Weinberg  and  Schweinler  (i),  in  1948,  were 
the  first  to  discuss  the  solution  for  a  reactor  driven 
by  an  oscillating  absorber.  Other  papers  (5,  5.)  have 
presented  results  for  neutron  wave  propagation  in  both 
multiplying  and  non -multi plying  media,  showing  the  dis- 
persion nature  of  the  process,  and  the  attenuation  and 
phase  shift  of  the  neutron  flux  disturbance. 

Further  evidence  of  the  spatial  effect  on  the  pro- 
pagation of  a  neutron  flux  disturbance  was  provided  by 
Badgley  (1)  and  Boynton  (8.).  In  their  work  with  random 
noise  measurements,  both  used  a  transfer  function  for  a 
two  point  reactor  model  to  represent  the  twin  slab  Uni- 
versity of  Florida  training  reactor.  Boynton  (8.)  was 
able  to  measure  the  transit  time  of  a  neutron  flux  dis- 
turbance as  it  travelled  between  the  two  slabs. 
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Statistical  correlation  techniques  were  used  as 
early  as  1946,  by  deHoffman  (9.),  to  measure  the  dimamic 
parameters  of  a  chain  reactor.  This  work  was  extended 
by  Albrecht  (10.)  and  Velea  (11)  to  include  the  reactivity, 
neutron  lifetime,  and  delayed  neutron  effects  in  the  the- 
oretical auto  correlation  function  of  the  neutron  density 
fluctuations. 

Moore  (i2,  13.)  showed  how  the  power  spectrum  of  the 
reactor  noise  is  related  to  the  reactor  transfer  function, 
where  the  power  spectrum  and  the  auto  correlation  function 
are  Fourier  transform  pairs.  Power  spectra  measurements  of 
the  reactor  noise  were  performed  by  Cohn  (14.)  and  Griffin 
(15 )  to  obtain  the  reactor  dynamic  parameters  as  expressed 
by  a  theoretical  transfer  function.  Griffin  v/as  able  to 
compare  the  noise  results  with  pile  oscillator  data. 

Certain  disadvantages  inherent  in  the  reactor  noise 
power  spectrum  and  auto  correlation  function  measurements, 
particularly  the  uncertainty  in  the  system  input,  the 
effects  of  external  instrument  noise,  and  the  lack  of 
phase  shift  information,  are  not  present  in  the  cross 
power  spectrum  and  the  cross  correlation  function.  Sev- 
eral investigators  (i,  ii,  17),  using  random  reactivity 
inputs  and  cross  correlating  the  input  and  output  signals 
of  the  nuclear  reactor,  have  measured  the  reactor  dynam- 
ic parameters  quite  easily,  even  with  input  signals  that 
were  small  compared  to  the  power  level  of  the  reactor. 


CHAPTER  II 
THEORY 

The  theory  presented  in  this  chapter  consists  of 
two  parts.  First,  the  general  spatially  dependent  trans- 
fer function  for  nuclear  systems  is  developed.  It  is 
then  adapted  for  several  special  cases  and  applications. 
It  is  written  for  a  point  source,  compared  with  the  common 
lumped  parameter  transfer  function,  and  converted  to  an 
alternate  form  that  converges  faster  for  large  systems. 
This  conversion  is  easily  made  for  a  non-multiplying 
medium,  but  for  a  multiplying  medium,  the  solution  of  a 
complex  transcendental  equation  is  required.  A  completely 
analytical  solution  in  terms  of  the  elementary  functions 
is  not  possible,  but  numerical  methods  may  be  used  to  ob- 
tain the  theoretical  transfer  function  for  specific  ap- 
plications. 

In  the  second  part,  the  statistical  theory  that  re- 
lates the  cross  correlation  and  cross  power  spectrum  of 
three  neutron  detectors  to  the  spatially  dependent  trans- 
fer function  is  developed. 

Transfer  Function 

The  theoretical  spatially  dependent  transfer  func- 
tion relates  the  thermal  neutron  density  fluctuations  at 
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the  point  r  in  a  nuclear  system  to  the  fluctuations  of  an 
artificial  source  of  fast  neutrons  having  a  known  spatial 
distribution. 

In  the  model  used,  the  nuclear  system,  which  is  iso- 
tropic and  homogeneous,  consists  of  a  single  region.  The 
driving  function  or  input  to  the  system  is  a  time  varying 
source  of  fast  neutrons.  The  fast  neutrons  experience  a 
slowing  down  phase,  represented  by  the  time  dependent  Fermi 
age  theory.  Upon  reaching  thermal  energy,  the  neutrons  en- 
ter a  diffusion  phase  at  constant  energy,  represented  by 
the  time  dependent  diffusion  theory.   If  the  medium  is  a 
multiplying  one,  additional  neutrons  are  created  by  fission, 
and  join  the  source  neutrons  in  the  slowing  down  process. 

It  is  assumed  that  the  neutrons  from  the  artificial 
source,  the  prompt  neutrons  from  fission,  and  the  delayed 
neutrons  have  the  same  energy  at  birth.  This  energy  is 
defined  as  the  zero  point  on  the  lethargy  scale. 

The  transfer  function  is  limited  by  the  same  restric- 
tions and  assumptions  inherent  in  the  Fermi  age  and  dif- 
fusion theories,  which  may  be  found  in  Reference  (19). 
and  by  the  mathematical  approximation  used  in  representing 
the  neutron  source  condition. 

Basic  eguatJons 

The  time  dependent  Fermi  age  equation  for  contin- 
uous slowing  down,  absorption,  and  leakage  iZj,    18.)  is 

1   t)(6(r,u,t)       'Dq(r,u,t)    r-     ^ 

V,   ■       =  -  >r -  Y^Jn)   0(r,u,t) 

v(u)  Ot  t)u  a    >-  '  ' 
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+  D(u)  \/^(t)iT,n,t)   +  S^(r,t)  6(u) 
5 
+  i;L:  e  (1  -  /})  (P^ir^t)    6(u)  +  Y.  AiCi(r,t)  5(u) 

i=l  (1) 

U  2  0 

Neutrons  starting  the  slowing  down  process  at  zero  leth- 
argy are  contributed  by  an  artificial  source  S^(r,t), 
by  prompt  fission,  VYs^e  H-fi)   <Pt(r,t),  and  by  the 
decay  of  the  neutron  precursors  A^C^(r,t).  The  slowing 
down  neutron  flux  is  0(r,u,t),  at  position  r,  at  lethargy 
u,  and  time  t,  in  units  of  neutron  /cm  sec  (unit  lethargy). 
The  neutron  velocity  is  v(u),  corresponding  to  lethargy  u, 
and  q(r,u,t)  is  the  slowing  down  density  at  r  and  t,  in 
units  of  (neutrons  passing  lethargy  u)  /cm  sec.  The 
macroscopic  absorption  cross  section,  L^(u),  and  the  dif- 
fusion coefficient,  D(u),  are  both  functions  of  the  neutron 
lethargy.  The  average  number  of  neutrons  produced  per  fis- 
sion is  symbolized  by  V,   the  thermal  fission  cross  section 
is  Ttf,   the  fast  fission  factor  is  €,  and  /)  is  the  frac- 
tion of  neutrons  produced  by  fission  that  appear  as  delayed 
neutrons.  The  thermal  flux,  0^(r,t),  is  not  a  function  of 
lethargy,  but  is  a  function  of  position  r  and  time  t,  with 
units  of  neutrons /cm^  sec.  The  six  neutron  precursors  are 
represented  by  C. (r,t),  i  =  1  to  5,  with  units  of  precur- 
sors/cm^, and  the  decay  constants  are  A.,  with  units  1/sec, 
The  Dirac  delta  function,  0(u),  is  used  to  insure  that  all 
the  neutron  sources  are  considered  mathematically  as  an 
initial  condition  in  lethargy. 

The  time  dependent  diffusion  equation  is  (1,  18.) 


,^t 


+  q(r,u^,t)  (2) 


with  the  symbols  previously  defined.  The  subscript  t  sig- 
nifies the  value  of  the  variable  at  thermal  energy.  -If 
desired,  the  model  used  here  could  be  changed  to  include 
an  artificial  thermal  neutron  source.  The  source  term 
would  be  added  to  the  diffusion  equation. 

The  equations  for  the  six  delayed  neutron  precursors 
are  (1,  18.) 


t)t 


1  =  1  to  6  (3) 


where 


i=l 

The  three  basic  equations  describing  the  neutron  be- 
havior are  coupled.  The  diffusion  equation  is  coupled  to 
the  slowing  down  equation  by  q(r,u^,t),  the  precursor  equa- 
tion is  coupled  to  the  diffusion  equation  by 
^^^' ^.^€0  .  (r,t ),  and  the  slowing  down  equation  is  coupled 
to  both  the  diffusion  equation  and  precursor  equation  by 
[t'Sfed./})  0t^r,t)  +  ^E^AiCi(r,t)]  5(u). 
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The  slowing  down  equation  may  be  written  in  terms 
of  q(r,u,t)  by  using  an  approximation  relating  <p{T,\x,t) 
and  q(r,u,t)  (18.) 

q(r,u,t)  =  £Et(u)  0(r,u,t)  (4) 

This  is  a  reasonable  assumption  if  ^g(u)  0(u)  is  slowly 
rarying  with  u,  and  Eg  (u)  =**  E^(u)  (18.), 

The  average  loss  in  the  logarithm  of  the  neutron 
energy  per  collision,  £,  has  units  of  (average  gain  in 
lethargy)  per  collision,  and  ^^^(u)  is  the  total  cross 
section  at  lethargy  u. 

The  relation  between  the  thermal  neutron  density  and 
the  thermal  flux  is  (18.) 

0^(r,t)  =  v^n^(r,t)  (5) 

Separation  of  time  varying  and  steady  state  compon- 
ents* The  time  varying  quantities  in  the  basic  equations 
may  be  separated  into  steady  state  components  and  time 
dependent  components.  The  latter  vary  positively  and  neg- 
atively about  the  steady  state  value.  Using  the  A  to 
represent  this  variation,  the  quantities  of  interest, 
after  substitution  of  Equations  (4  and  5)  in  the  first 
three  equations,  are 

nt(r,t)  =  \o^^^   +  An^(r,t) 

q(r,u,t)  =  qo(r,u)  ^-  Aq(r,u,t) 

C.  (r,t)  =  C.^(r)  ^   Ac^(r,t) 

S«(r,t)  =  S,^(r)  +  As  fr,t)  (6) 
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Since  the  basic  equations  are  linear  with  respect 
to  the  time  varisible,  they  may  be  separated  into  the 
steady  state  and  the  time  varying  equations. 

For  the  steady  state  components,  the  slowing  down 
equation  is 

t)q  {r,u)    E  (u)  q  (r,u)     D(u)     « 

+  S^^(r)  6(u)  +  VL^ea  -  J3)   v^n^^Cr)  5(u) 

6 
+  Y.   AiC.^(r)  6(u)  (7) 

i=l 

the  diffusion  equation  is 

0  =  -  JlatVto^^'  +  VtV%o^^^  ^  '^fo^^'^t'   (8) 
and  the  precursor  equations  are 

0   =  A^^f^-t^to^^^  -  ^i^io^-^ 

i   =  1  to  6  (9) 

The  time  varying  components  satisfy  the  following 
equations.  The  slowing  down  equation  is 

1     ^Acr(r,u,t)     'OAq(r,u,t) 


^Et^^J  ^^^^   ^t  t)u 

E„(u)    .  D{u)     «  . 

-  -r^ Aq(r,u,t)  +  ~-=r-— -  U  Aq(r,u,t) 

£5:,(u)  £Z,(u)  V 

+  As^(r,t)  6(u)  +  I'T.^ea   -  /))  v^An^Cr^t)  6(u) 
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+  ^AiAc,(r,t)  6(u)  tlO) 

i=l 


the  diffusion  equation  is 

t)t  at  t    t 

+  Aq(r,ut,t)  (11) 

and  the  precursor  equations  are 

^^^i^"'^^.  =  /),l^E,ev,  An,(r,t)  -  A,Ac,(r,t) 

i   =  1  to  5   (12) 


Solutions  of  time  dependent  equations 

Equations  (10,  11,  and  12)  are  now  solved  for  the 
transfer  function,  ^ghich  is  defined  as  the  variation  or 
perturbation  in  the  thermal  neutron  density,  An^(r,(J), 
at  position  r,  caused  by  a  variation  in  the  artificial 

source,  A^^^^'^^  6(u). 

Laolace  transform  in  time.  A  Laplace  transforma- 
tion of  Equations  (10,  11,  and  12)  is  now  performed,  which 
transforms  the  equations  from  the  time  domain  to  the  fre- 
quency domain  and  removes  the  time  derivatives  from  the 

equations.  The  time  varying  components  are  restricted  to 
bounded  functions,  which  may  be  composed  of  sinusoidal 

components.  Thus  the  Laplace  transform  converges  with  the 

real  part  of  the  Laplace  variable,  s  =  CT  +  jCJ,  equal  to 

»erto,  or  s  =  j  (J  iW . 
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Using  the  initial  condition  that  all  the  time  de- 
pendent variables  are  zero  at  time  equal  zero,  which 
eliminates  the  transient  solutions.  Equations  (10,  11, 
and  12)  become  respectively 

jUAq(r,u,(J)  _    t)Aq(r,u,  U)       S^(u)  Aq(r,u,(J) 
|Lt(ti)  v(u)   "  "   t)u         '      ^I^t^'^^ 

V^  Aq(r,u,U)  +  As  (r,U)  6  U) 


£Lt(u) 

+  y  ZjCd  -  /))  v^  Ant(r,U)  5(u) 

6 
+  ^  A.  Ac^(rA>0  (5(u)  (10a) 

i=l 

j(jAnt(r,(J)  =  -  Z^^.v^  Ant(r,U)  +  D^v^^^  ^  n^  (r,U) 

f  AQ(r,ut,U)  (11a) 

jU  ACi(r,U)  =  /JiZ^^f^v^  /^n^ir.U)    -  A^  ACi(r,U) 

i   =  1  to  6   (12a) 

The  independent  time  variable,  t,  is  replaced  by  an 
independent  frequency  variable,  CJ,  in  all  the  terms. 

Spatial  expansion.  The  slowing  down  density,  the 
thermal  neutron  density,  the  precursor  density,  and  the 
artificial  source  are  now  expanded  in  terms  of  complete, 
orthogonal  functions  that  satisfy  the  Helmholtz  equation 
ill) 


V  ^^^'V     =  -  S  '^^'"'^P^ 


(13) 
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over  the  volume  of  the  system.  The  function  H(r,B  )  may- 
be sines  and  cosines  for  Cartesian  coordinates,  and  for 
cylindrical  coordinates  the  function  may  be  Bessel  func- 
tions and  sines  and  cosines.  The  requirement  that  the 
slowing  down  density,  thermal  neutron  density,  and  pre- 
cursor concentration  be  equal  to  zero  at  the  extrapolated 

boundaries  of  the  system  is  satisfied  by  H(r,3  ).  The 

2 
eigenvalues,  B  ,  which  correspond  to  the  familiar  buckling 

in  reactor  analysis  terminology,  may  be  composed  of  one 

component  for  each  space  dimension;  for  example 

in  general,  or 


■n,ni^     rn,7ri^  rn^TTi^ 


■; '  ra  -m  i 


2c  J 
n- ,  n-,  n  all  odd   (14a) 

for  a  rectangular  parallelpiped. 

Performing  the  expansions,  the  four  freqfuency  depen- 
dent variables  become 

Aq(r,u,U)  =  ^qp(u,(J)  H(r,Bp) 
P 

[\n^{r,U)       =  ^np(U)  H(r,Bp) 
P 

Ac,(r,U)   =  Y.%^^^  ^"^^""'V 
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diS^ir.U)      =  ^3p(CJ)  H(r,Bp)         (15) 
P 
The  bar  symbolizes  the  expansion  coefficients. 

To  solve  Equations  (10a,  11a,  and  12a)  for  the  trans- 
fer function,  it  is  necessary  that  the  spatial  distribu- 
tion of  the  artificial  neutron  source  be  known.  Therefore, 
^  may  be  found  by  using  the  orthogonality  property  of 
H(r,Bp)  (19). 

j/\\iT,U)   H(r,Bp)  d(Vol)  =  J  ^Sp(CJ)  H(r,  B^) 
Vol  Vol  p 

(•)  H(r,Bp)  d(Vol) 
=  ^p^^^  \  (16) 

since  (19) 

jH(r,Bp)  H(r,Bp)  d(Vol)  =  Ap 


Vol 


B   =  B 
P      P 


Bp  #  Bp(17) 


A  is  a  constant  dependent  on  p  only. 

Substituting  Equation  (15)  into  Equations  (10a,  11a, 
And  12a),  and  using  Equation  (13)  where  applicable,  the 
•slowing  down  equation  becomes 
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r-,       r  _      f    J  (J      ^a(^) 
7H(r,B  )   -q  (u,U)  -T-?: +  "14= 


D(u)  ^l\       tiq^Cu^U) 


^-  I'LfCd  -  /))  v^flp((j)  4-  ^p(U)|  (5(u)l    =  0  (18) 

the  diffusion  equation  becomes 

^H(r,Bp)   [.np(U)    {iU  *    L^^t  *  Vt^^} 

P 

+   qp(ut,CJ)]      =     0      (19) 
and  the  precursor  equations  become 

^H(r,Bp)  [-Cip(U)   {jU  +  A.}   +   P.VL^er^n^iU)]     =    0 
P 

i      =     1  to  6      (20) 

3y  use  of  Equation  (13),  the  spatial  divergence  term 
has  been  replaced  by  -3  H(r,BL)  from  the  Helmholtz  equa- 
tion. Equation  (18)  is  the  only  one  still  containing  a 
derivative. 

Since  the  sums  in  Equations  (18,  19,  and  20)  must 
equal  zero  independently  of  r,  each  coefficient  of 
H(r,B  )  must  be  equal  to  zero  (19.).  Thus,  the  terms  en- 
closed by  the  square  brackets  can  be  set  equal  to  zero 
individually. 

Slowing  down  solution.  The  slowing  down  equation  has 
now  been  reduced  (from  Equation  (18))  to 
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qp(u,CJ) 


jCJ 


E^ 


D(u)  Bj 


£^(u)J 


^Et(^)  ^(^)    £^t^^^ 


,(U)  +  y  Ef  Cd  -  P)  v^Hp(U) 


■f-    E'^i^ip(cj) 


1=1 


6(u)  (21) 


viiiere  the  bracket  multiplied  by  0  (u )  is  the  initial  con- 
dition on  q  (u,CJ). 

The  solution  to  this  equation,  subject  to  the  boundary 
conditions,  is  easily  found  and  is 


qp(u,U)  =  qp(0,L;)  exp 


/     du 
-J(J   -TF 


-I 


i^E  (u)  du   „  ;^D(u)  du 


-B^ 


^L,M     pJ£Et(u)J 


/ 


(22) 


The  first  integral  in  the  exponential  function  can  be 
identified  as  the  average  time  required  to  slow  down  from 
jsero  lethargy  to  lethargy  u.  This  time  is  defined  as 
Jirr'iyi),   with  the  slowing  down  lifetime  to  thermal  lethargy 
defined  as 


i^K^  =   ^r 


(23) 


The  exponential  function  raised  to  the  power  repre- 
sented by  the  second  integral  is  known  in  reactor  analysis 
as  the  resonance  escape  probability  P(u)  (18.). 


16 
The  third  integral  is  defined  in  reactor  analysis  as 
the  Fermi  age.  Tin),    (li).  The  Fermi  age  to  thermal 
lethargy  is  defined  as 

'^K^  =  '^  (23a) 

Thus,  Equation  (22)  becomes 

qp(u,(J)  =  qp(0,U)  P(u)  ezp[.jU/^(u)  -  Tin)   b2](24) 

The  slowing  down  density  at  zero  lethargy,  qp(0,(J),  may 
be  found  from  the  initial  condition 


qp(o,(j) 


%(U)  +  vL^ea  ~  J3)  v^n^iu) 


i=l 


Pyecurspr  tprmp.  A  relation  between  the  precursors 
and  the  thermal  neutron  density  is  provided  by  Equation 
(20) 

=   A^^f^v.n,(U) 

j(j  +  A. 


^ip(U)   = 


(26) 


The  number  of  unknown  variables  in  the  basic  set  of  equa- 
tions may  be  reduced  by  substituting  this  relation  into 
Equation  (25),  obtaining 


qp(o,c.;)   =   Sp(cj)  +  i^^fCv^ 


(27) 
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For  convenience.  Equation  (27)  may  be  simplified  by 
defining 


5 


J3.X. 


Slowing  dovm  density  terms.  Substitution  of  Equation 
(27)  into  Equation  (24),  then  Equation  (24)  into  Equation 
(19)  further  reduces  the  basic  equations  by  eliminating  the 
islowing  down  density,  q  (u,CJ),  resulting  in 

+  VLf^J3^v^R^iU)]    P(ut)  exp[-jU/^  -  TB^]   =  0   (29) 

Solving  for  n    (CJ),    and  using  the  following  defini- 
tions   (18.) 

2  \ 

L  =^^ 


^at^t 

yZpP(uj  e 


k       =    1>— 1 (30) 

00  y 

^at 


for  the  square  of  the  thermal  diffusion  length,  thermal 
lifetime  in  an  infinite  medium,  and  multiplication  factor 
in  an  infinite  medium,  respectively,  gives 

-  „  „  _  %<^^  PKI  4 

n-,((J)  - 


(31) 
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General  transfer  function 

An  expression  for  the  thermal  neutron  density  may  be 
constructed  using  Equations  (15  and  31). 

(32) 
S  (CJ)  is  found  using  Equation  (16)  and  fj^   is  defined  by 
Equation  (28). 

This  is  the  general  transfer  function,  showing  the 
perturbation  in  the  neutron  density  at  location  r  caused 
by  a  source  A  3^(r,Cu/). 

General  point  source  transfer  function 

The  transfer  function  can  be  obtained  in  a  more  ex- 
plicit form  if  the  artificial  neutron  source  is  a  point 
source  located  at  r^.  'ilien,  from  Equation  (16) 

^p<^^  Ap  =  J  As^(r,CJ)  6(r  -  r^)  H(r,Bp)  d(Vol) 
Vol 

=  As^(ro,(J)  H(ro,Bp)  (33) 

The  correct  form  of  the  Dirac  delta  function  corresponding 
to  the  geometry  of  interest  must  be  used. 
Therefore,  Equation  (32)  becomes 
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=     S(r,r^,U) 


=  ^K^-fS: 


H(r„.B^)  H(r,B„) 


*   "Tb^)-  k„/3J    (34) 


Thus,   the  perturbation  of  the  thermal  neutron  density  at   r, 
caused  by  a  fast  neutron  source  at  r^,    is  equal  to  an  in- 
finite sum.     Each  term  of  the  sum  is  composed  of  the  pro- 
duct of  the  orthogonal   eigenf unctions,   divided  by  the 
characteristic  equation  for  the  Fermi   age,   diffusion 
theory  model  used  in  this   study.     The  characteristic  equa- 
tion corresponds  to  the  dispersion  function  of  classical 
physics.     Vftien  set  equal  to   zero,   this   function  is   called 
the  dispersion  law   (42).     Equation    (34)   is  also  the  Green's 
function  for  this  model,   and  could  be  used  as  a  kernel   in 
integrating  over  a  source  occupying  a  finite  volume    (27). 
Equation    (34)   can  be  written  in  another  form,  using 
the  more  easily  measured  quantities 

•='  (35) 


1  +  l^bI 


t     exp(-TBfj  ^      ^ 

®^^  1  +   L^BJ 


which  are  the  thermal  neutron  lifetime  and  multiplication 
factor  for  a  finite  medium. 
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Then 

G(r,r^,(J)     =     [P(ut)  J?g  exp(-rBf)] 


■'Ik 


H(r^,B^)  H(r,Bp) 
2/qZ    _    q2 


i.iii!ijjii.jcj£ 


1  +  L^B^ 


exp 


[j(Ji^ 


r(B^  -  b2)]  .  k^,,/)  1 


(37) 


There  are  several  interesting  features  about  the 
transfer  function  as  vrritten  in  Equations  (34  and  37). 
The  denominator,  ^^Axich  is  the  characteristic  equation  of 
the  system,  is  readily  available  for  stability  and  control 
studies.  The  transcendental  equation  can  be  solved  numeri- 
cally for  its  roots,  or  an  approximate  form  can  be  used 
(see  Appendix  C). 

The  approach  used  above  to  find  the  expansion  coef- 
ficients for  the  source  may  still  be  used  even  if  the 
source  is  located  on  one  face  of  the  assembly.  The  eigen- 
function,  H(r,B  ),  would  then  be  defined  over  the  actual 
medium  and  its  mirror  image,  putting  the  source  at  the 
center  of  the  mathematical  problem. 

The  quantity  defined  by  Equation  (28),  JJc  is  simpli- 
fied at  the  two  extremes  of  the  frequency  range.  As  the 
frequency  goes  to  zero,  /J ^  approaches  one.  As  the  fre- 
quency goes  to  infinity,  J3^  approaches  1  -  /j.  Only  in 
the  frequency  range  0.001  cps  to  0.5  cps  is  JJ^  changing 
significantly.  Thus,  /}q  may  often  be  replaced  by  its  lim- 
iting value,  depending  on  which  frequencies  are  of  interest, 
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Both  forms  of  the  transfer  function.  Equations  (34 
and  37),  have  the  disadvantage  that  the  convergence  of  the 
series  is  slow,  particularly  if  the  size  of  the  system  is 
greater  than  the  diffusion  length  and/or  the  Fermi  age, 
Since  one  of  the  assumptions  of  the  Fermi  age  and   diffusion 
theories  {JIS.)  is  that  the  system  is  large  compared  to  the 
characteristic  nuclear  dimensions,  many  terms  must  be  in- 
cluded for  a  numerical  analysis  of  Equation  (34)  or  Equa- 
tion (37).  An  alternate  form  of  the  transfer  function  that 
converges  faster  as  the  size  of  the  system  increases  is  de- 
veloped after  the  next  section. 

Comparison  with  the  lumped  parameter  transfer  function 

It  is  interesting  to  note  that  Equation  (34)  can  be 
readily  converted  to  the  lumped  parameter  transfer  func- 
tion commonly  used  for  nuclear  reactors. 

First,  let  all  the  quantities  associated  with  the 
Fermi  age  model  go  to  zero,  that  is,  T  to  0,  u.  to  0,  and 
£^   to  0.  Second,  all  terms  in  the  series  beyond  p=l  may 
be  ignored,  or  better  yet,  it  may  be  assumed  that  the  in- 
tensity of  the  artificial  source  is  given  by 

A3a(r,(J)  =  A3(CJ)  H(r,B^)  (38) 

which  is  the  fundamental  mode  of  the  system.  Then  ^  is 
zero  for  all  p  greater  than  one,  and 
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As   (U)  '1+  B^L^^  jUH^  -  k^/}^ 

H(r  B^) 


(39) 


If  the  thermal  neutron  density  is  measured  at  the 
center  of  the  system,   where  H(0,B.)=1,    then  Equation    (39) 
is  identical  with  the   lumped  parameter  transfer  function 
(5.)   for  a  subcritical   system;    otherwise,   they  differ  by 
a  constant.     As  k^ff  goes  to  one,   the  correct   form  for  a 
critical  reactor  is  obtained 

An(r,(J)  H(r,B-) 

^  1- ; _ (40) 


As  iU)  ,  ,      1 


AjA 


Alternate  form  for  the  t>oint  source  transfer  function 

The  general  transfer  function  may  be  converted  to 
another  form  that  converges  much  faster  for  large  systems. 
The  procedure  followed  is  known  as  the  inversion  theory  of 
series  (2£),  and  as  Poisson's  summation  formula  (4,).  In 
general,  the  method  may  be  applied  in  three  dimensions  in 
any  coordinate  system.  However,  for  the  specific  approach 
described  below,  it  is  necessary  that  the  eigenvalue  of 
Ecjuation  (13)  be  periodic. 

gpecialized  problem.  In  order  to  show  the  inversion 
in  some  detail,  the  geometry  of  the  medium  is  specified  as 
cylindrical,  the  same  geometry  as  that  used  in  the  experi- 
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ments  (see  Chapter  IV).  Since  the  axial  eigenvalue  is  the 
only  one  that  is  periodic  in  cylindrical  coordinates,  the 
inversion  is  performed  in  one  dimension  only. 

The  artificial  neutron  source  is  located  on  the 
centerline  of  the  experimental  assembly  and  approximated 
by  a  point  source.  The  origin  of  the  coordinate  system 
is  located  at  the  source. 

For  the  cylindrical  geometry  described  in  Chapter  IV, 
the  eigenfunction  is 

mTTz 


(41) 


with 


3^ 

B^ 

"^' 

2 

+ 

'mjT' 

_  R  _ 

_2c  _ 

m  =  1  to  OD  ,  odd  only 

n  =  1  to  00  (42) 


The   variable  P  is  the  radial  coordinate  divided  by  the 
maximum  radius,  R.  The  n^^  zero  of  the  J^  Bessel  function 
is  V^,   and  c  is  the  distance  from  the  source  to  the  end  of 
the  assembly  in  the  positive  z  direction,  even  if  the  source 
is  on  one  face  of  the  cylinder.  Under  the  assumptions 
listed  above,  the  solution  is  valid  only  for  0 -=  z  <  c  if 
the  source  is  not  in  the  center  of  the  medium. 

Using  Equation  (16)  to  find  iL  and  ^  ,  Reference  (21) 
gives  for  cylindrical  geometry 
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\> 


=     R^ 


I^o^VnP^  Jo^KP' 


'TT: 


LVol 


cos  

2c  2c 


(•)   2TCpdPdz 


For  the  point   source  at   r  =  0, 


(43) 


Apsrp(u)  =   r 


Ta         <5(p) 


6(z)  J^(^'^p) 


.Vol 


(•)   cos  27tPdPdz 


2c 


therefore. 


=      27rR2A3((J) 


^^(CJ) 


2  AS(U) 


and  Equation    (32)  becomes 

Ant(r,CJ)      _      2P(u^)  jg^ 


(44) 


(45) 


G(r,o,U)    = 


(•) 


As(CJ) 


c 

mTTz 


^^ 


Z^  Z^     J?  ( 1/  )     f(i  ♦  lV  +  jU/  )  sxpUUIt 

n=l  m=l   1    "   L       P      =         ' 

tt 

where  sin  (^^^)  has  been  included  in  the  numerator  so  that 
m  may  take  all  values,  from  one  to  infinity. 
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2 
Inversion  procedure.  The  first  step  is  to  replace  B^^^ 

2 
in  Equation  (46)  by  a  continuous  variable,  ^  .  To  com- 
pensate for  this  replacement.  Equation  (46)  nnist  now  be 
multiplied  by  the  periodic  delta  function 


y   6(^-2IL)  (47) 


The  sum  over  the  index  m  is  replaced  by  one  half  of  the 
integral  over  c    from  minus  infinity  to  plus  infinity. 
Equation  (46)  becomes 


P(u.)0     -r-  JJV^P) 
G(r,o,(J)     =     ^-^    )     -| J^^ 

n 


Llfa\ 

d^    cos£  z  sin^  £°    E    ^  ^£    "  2^  ' 


(•) 

-<D 


r        S         S  b      ^_    ^   :>        ^c 

I      cr=-oD 

J        [(1  +    L^l    +l2£2^   J(J4)exp(j(J£^ 

The  periodic  delta  function  may  noi^  be  replaced  by  its 
Fourier  expansion.  It  is  shown  in  Appendix  A  that 

Z<5(£    -f?'     =     f   fexp(j4mo|)  (49) 

qr-00  m=:-oo 

Assuming  that  the  summat^'on  over  m  and  the  integration 
commute,   their  order  may  be  interchanged.     The  integral   is 
seen  to  have  the  form  of  a  Fourier  transform  in  the     £" 
variable. 
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Since  Equation    (48)   is  an  even  function  of  the  ^vari- 
able,  the  cosine  functions  may  be  replaced  by  expCj  6),    or 
vice  versa.      Doing  this,    the  summation 

00 

^^^^     "       y    cos£  z  sin^^c  cos   4mc£  (50) 

m=-oo 

may  be  reduced  to  a   sinq^ler  form  as   shown   in  Appendix  B, 

The  result   is 

00 

f(^)    =      7    ^      (.1)°^  cos£(a  +  2cm)  (51) 

m=-oo 

00 

-      i    E      ^"-^^"^  exp[j^(z  +  2cm)]  (52) 

m^-OD 

Equation    (48)  becomes 

00 


) 


P{u.)   i'.V"   ■'o^*'i 


a=l     ^"'- 


Z°°  „  r* 9xp[jU(z  +  2cm)]  d^ 


+  rB?  +  rf)  -^c.fic]  '" 


Thus,  the  original  cosine  eigenf unctions  have  been  replaced 
by  positive  and  negative  "sources"  located  periodically 
along  the  z  axis. 

The  integral  in  Equation  (53)  may  be  evaluated  by  the 
residue  theorem  of  complex  variables  (19.),  as  applied  to 
real  integrals.  The  main  steps  in  performing  the  integration 
are: 
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1.  The  real  variable  £  is  replaced  by  the  conq^lex 
variable  Q  . 

2.  The  poles  of  the  integrand  are  found,  corresponding 
in  this  case  to  the  roots  of  the  denominator  of  Equation  (53) 

DP  =  (1  +  l2b2  -I-  l2^2+  iUf^)   exp(JU/^  +  T 3^ 

*rt,^)  -^^fic     (54) 

3.  The  residues  are  evaluated  at  the  poles  in  the 
upper  half  plane  when  z  f  2cm  is  positive  and  at  those  in 
the  lower  half  plane  when  z  ♦•  2cm  is  negative.   The  variable 
z  is  restricted  to  values  between  o  and  c. 

4.  The  value  of  the  integral  is  equal  to  the  sum  of 
the  residues  multiplied  by  2TCj. 

The  roots  of  the  denominator  cannot  be  found  in  terms 
of  elementary  functions,  since  the  denominator  is  a  trans- 
cendental equation  having  an  infinite  number  of  complex 
roots.  Representing  these  roots  by 

^i  =  ±JRi  =  ±j{A.  +  jB^) 

i  =   -oo  to  CD   (55) 
allows  a  symbolic  solution  of  the  integral,  with  the  small- 
est value  of  R^  corresponding  to  i  =0.  Appendix  C  pre- 
sents a  discussion  of  the  characteristic  equation,  its 
roots,  and  how  thev  may  be  found. 

The  residue  contributed  by  the  roots  of  DF  are 

exp(-JtJA-  -    Tb^  +  r^f)   exp[-R.  (z  +  2cm)] 
Res.     =     r-K* 5—5 *   A  '  1 vob; 

2jI?i[L^   +    7(1+   LX   -  LRi+   J<*^4)] 
where  L'Hospital's   rule  has  been  used.      It   is  possible  for 
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a  root  of  DF  to  lie  on  the  real  axis,    if  k^ff  ^  1.     For  this 
case,   the  residue  is  ecpial  to  one  half  the  value  obtained 
using  Equation    (56). 

gumming  the  residues,   the  transfer  function  becomes 

°°    Jo^^nP^ 


G(r,o,(J)     =     P(ut)j?s     exp(-jUjP^)    ^    jl^y""]' 

n=l 

i=?-oo       ■• 

+  Y     (-1)°^  {gxp[-R.  (z  +  2cm)]  +  exp  [h.  (z  -  2cm)]}     (57) 
m=l 

Thus,  the  transfer  function  for  a  point  source  in  a 
cylindrical  medium,  with  a  detector  located  at  r,  is  equal 
to  the  product  of  the  sum  over  n,  the  radial  coordinate 
index;  the  sum  over  i,  the  roots  of  the  characteristic 
equation  (Equation  (54))  in  the  upper  half  plane;  and  the 
sum  over  m,  the  position  index  on  the  z  axis. 

The  transfer  function  (Equation  (57))  for  systems  that 
are  large  compared  to  the  characteristic  nuclear  dimensions 
(T  and  L^)  will  ordinarily  be  adequately  represented  by 
just  a  few  of  the  possible  number  of  terms. 

Alternate  form  for  the  non 'multiplying  medium  transfer 
function 

For  the  case  where  k^^  is  equal  to  zero.  Equation  (45) 
may  be  easily  inverted,  producing  a  result  written  in  terms 
of  the  elementary  functions.  The  problem  will  be  restricted 
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to  a  cylindrical  geometry  containing  a  point  source,  as 
discussed  in  the  previous  section. 

Inversion  Procedure.  The  procedure  followed  in  the 
previous  section  is  applicable  to  the  case  for  which  k^^ 
is  equal  to  zero.  Thus  Equation  (53)  may  be  re-written 


n=a 


00 


(58) 


The  roots  of  the  denominator  are  easily  found,   and  are 


Li 


(59) 
(50) 


where 


O^n     = 


[(1  4-   lV)^-!-    U^IY  ^  ^  "•  ^^^n 


2V 


(61) 


and 


r  = 


[d  +   L^B^)^  +   CJ^ig]      -    (1  +  L^B^) 


2V 


(52) 


Integrate  ng  around  the  upper  half  plane  when  z   +  2cm 
is  positive,   and  the  lower  half  plane  when   z   +  2cm  is 
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negative,  the  transfer  function  becomes 


G(r,o,(J)  =   V^   e^ 


Tf^i^^^T^-  h^ 


n=l  •- 

+  Y.     ^-^^''  exp[-(z  +  2cm)  iOCn   +  J  Tn^] 
m=l       ^  ^ 

♦  exp[(+z  -  2cm)  (a^  +  J  7n^]\ 


(63) 


Equation  (53)  applies  for  a  non -multi plying  medium  contain- 
ing a  fast  source  located  on  the  centerline  of  the  assembly 
and  for  0  «=  z  «:  c. 

Appendix  D  contains  a  description  of  the  con^suter  code 
used  to  evaluate  Equations  (57  and  63),  and  Chapter  V  con- 
tains the  theoretical  results. 


Statistical  Theory 

The  statistical  correlation  and  power  spectrum  equa- 
tions used  in  obtaining  and  processing  the  experimental  data 
are  developed  in  this  section. 

If  a  nuclear  system,  represented  by  its  spatially  de- 
pendent transfer  function,  is  driven  by  a  random  input  (the 
fast  neutron  source),  the  output  is  also  a  random  function. 
The  statistical  relations  between  the  output  and  the  input, 
or  between  several  outputs,  may  be  expressed  in  terms  of 
the  system's  transfer  function. 
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In  this  study,  the  cross  power  spectrum  is  calculated 
from  the  cross  correlation  of  two  outputs,  -which  are  the 
neutron  density  variations  at  two  different  locations  in  the 
medium.  Repeating  the  experiment  with  one  output  relocated 
and  using  the  equations  developed  in  this  section,  the  ex- 
perimental transfer  function  of  the  medium  can  be  obtained. 

Definitions 

There  are  several  definitions  that  are  basic  to  the 
statistical  analysis  of  random  functions.  One  that  is  very 
in?50rtant  is  the  cross  correlation  function  (22) 

-T 

It  is  seen  that  this  is  just  the  average  over  all  time  t  of 
the  product  of  two  functions,  with  one  function  delayed 
relative  to  the  other  by  the  time  T.  The  cross  correlation 
function  is  frecyuently  written 

^i2^r)  =  lYrrrTjrt-rrT  (ss) 

to  signify  the  averaging  process.  For  a  stationary  pro- 
cess, the  time  average  may  be  replaced  by  an  ensemble 
average  if  desired  (22) . 

The  Fourier  transform  of  the  cross  correlation  func- 
tion is  (22.) 


$i2(^>  =     nr     I    ^12^'^^  exp(.j(^r)  df 


(56) 
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and  is  defined  as  the  cross  power  spectrum  density. 

An  important  relationship  in  system  analysis  is  the 
convolution  integral,  which  can  be  used  to  relate  the  sys- 
tem characteristics  and  the  input  to  the  output  of  the  sys- 
tem (22) 

CO 

f^(t)  =  r  hiV)  fiCt  -  V)  dv        (67) 

-00 

where  f  (t )  and  f.(t)  are  the  input  and  output  functions, 

o  1 

respectively,  and  h(L/)  is  the  impulse  response  function  of 

the  system  iZZ) • 

The  impulse  response  function  and  the  system  frequency 
function,  H((J),  of  the  system  are  related  by  the  Fourier 
transform,  and  they  form  a  Fourier  transform  pair  (22) 


h(t)   =  ^  /  ^^(^>  exp(jCJt)  dU 
-co 

00 

H((J)  =   j    hit)   exp(-jCu't)  dt 


(68) 


For  the  class  of  functions  considered  here,  the  system 
frequency  function  and  the  transfer  function,  defined  by 
the  Laplace  transform,  are  essentially  the  same  thing, 
since  h(t)  is  equal  to  zero  for  times  less  than  zero,  and 
the  real  part  of  the  Laplace  variable  is  zero. 

Equation  (57)  may  be  used  to  show  the  relationship 
between  the  cross  correlation  of  the  input  and  output  of 
a  system  and  its  impulse  response  function.  If  fj^it)   and 
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f2(t)  in  Equation  (64)  represent  the  input  and  output  sig- 
nals, fj^(t)  and  ^^(t),  respectively,  then  substitution  of 
Equation  (67)  for  fo(t)  gives 

T  CD 

4^io^'^^  =  T^i^'oD  h     fdt  fi(t)  J  hiV)   fi(t 

+  r  -  V)  dV  (69) 


Assuming  that  the  two  integrations  commute,    their  order  may- 
be interchanged.     Realizing  that  Equation    (54)  also  defines 
an  auto  correlation  function  if  f2(t)  and  ^2^^^  ^^®  ^^® 
same.   Equation    (69)  becomes 


4^io^'^^      =      J    h(^')  Cj)^.(r    -    V)  dV 


(70) 


This  equation  has  the  same  form  as  Equation  (67),  but  re- 
lates the  statistical  correlation  functions  rather  than  the 
actual  input  and  output. 

The  Fourier  transformation  of  Equation  (70)  may  be 
performed  by  multiplying  both  sides  of  the  equation  by 

exp(-jCJt)  dT/27r  and  integrating  over  all  time  T.  Using 
Equation  (68)  for  H(CJ)  gives 

(^^^^iU)    =   H((J)  $i.((J)  (71) 

Equation  (71)  relates  the  system  frequency  function 
and  the  auto  power  spectrum  of  the  input  to  the  cross 
power  spectrum  of  the  input  and  output  of  the  system  (J2.). 
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Cross  power  spectrum  between  two  outputs 

Referring  to  Figure  1,  the  experimental  setup  is  re- 
presented by  an  input  function,  3(t);  three  spatially  de- 
pendent transfer  functions  between  the  source  and  detectors 
1,  2,  and  3;  three  measuring  system  transfer  functions;  and 
a  black  box  that  performs  the  cross  correlation  of  the  in- 
coming functions.  See  Chapter  IV  for  a  description  of  the 
experimental  equipment. 

Functions  f^,  f**  and  f-  represent  the  neutron  density 
variation  at  each  of  the  detectors  ^en  the  variation  is 
caused  by  the  source,  while  f^,  f^,  and  fg  are  uncorrelated 
variations  in  the  neutron  density  at  the  detectors,  or  noise 
as  far  as  the  experiment  is  concerned.  Functions  fy,  fg, 
and  fg  are  instrument  noise  added  to  the  signed  as  it  passes 
through  the  measuring  system.  Functions  f^o,  f^x*  ^^d  f22 
are  the  resultant  signals  that  are  cross  correlated.  From 
this  systepa,  the  experimental  transfer  function  is  obtained 

Ho«D.(U)  =   T^  ^  ^  ^ (72) 

Using  Equation  (57),  the  neutron  density  variation  at 
each  detector  is 

foi(t)  =  f^{.t)  ^   f4{t) 

fh^L^^^  ^^'^   -  u)  du  +  f4(t)       (73) 
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1' 


Source,  S(t) 


ml 


m2 


m3 


^10    ^11 


\f  ^ 


^10    ^12 
0DlD3^r) 


Figure  1.  Definition  of  functions  for  statistical  theory 
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CD 

[hgC  v)   3(t    -  V    )   dv  -H   fgCt)  (74) 

-00 

fD3    (t)    =      f3(t)  ^  fgCt) 

00 

fhgCw)  S(t  -  w)  dw  ^-  fg(t)  (75) 

-00 

The  signals  -which  are  cross  correlated,    f-j^Q,    f^.!'   ^"^ 
f,«,    are  affected  by  the  impulse  response  function  of   the 
measuring  system  and  are 

CD 

f;^0(t)     =     f7(t)      ^     j^^(^)   dx  foi(t   -  x) 

-00 

OD  - 

=     f7(t)     +       f^al^^^   ^  I  ^4^^   -  ^^ 

-OD 

00  -. 

+         r  hjL(u)   du  3(t    -  X   -  u)  I    (75) 
-oo 

00 

f3Li(t)      =     f8(t)     i-    j  h^2(y)   dy[f5(t   -  y) 
-co 

00 

+        j  h2(v)   3(t   -  y  -  v)   dv        (77) 
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00 


fl2(t)    =   f9(t)    +     J  h^3(2)  dJfgd 

mfVS  *• 


-CD 

00 


+      J  hgCw)   3{t   -  z   -  w)  dw  I     (78) 


-oo 

Thus,    the  functions  that  are  cross  correlated  are  related 
to  the  system  input  through  two  impulse  response  functions, 
one  for  the  nuclear  system  and  one  for  the  measuring  system, 
and  to  the  noise  present   in  the  two  systems. 

Equation    (54)  may  be  used  to  obtain  the  two  cross 
correlation  functions 


*D1D2    (T)     =     T^^oD    k  J    Wt)   fn(t    +  T)  dt      (79) 


,T 
-T 
and 


,T 

:t 


^DlD3<r)      =    T^i'oo   h  J     fio(t)   fl2(t    +    T)  dt    (80) 


Sttbstituting  Equations    (75  and  77)   in  Equation   (79) 
results  In 

T 


0(r )    =    ^^^    J»-  r  dt  If  (t) 

-T 

OD  OD  .  - 

+     Jhjjj^(x)  dx  [f4(t   -  x)    +    J    hi(u)   du  S(t    -  X  -  u)   J 
-00  '  -00 

[CD 
fg(t   +    T)    +    J  h^2(y)  dy     f5(t  +    T    -  y) 
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00 


+    r  hgCv)  dv  S(t    +    r   -  y  -  v)     j 


(81) 

•00 


Performing  the  multiplication,  separating  terms,  and  using 
the  definition  of  the  cross  correlation  function  (Equation 
(55)),   gives 


0DlD2(r)      =     fy^t)    fglt    ..    T) 

00 

+  f   \lU)  dX    f^Ct')     fg(t'+       f      +     X) 

-00 
OD 

+    /^ms^^^  ^y  f7(t)  f5(t  4-  r  -  y) 

-00 

OO  00 

+        rh^^(x)  dx      J    h^(u)  du  3(t')   f8(t'+    f    +    X  +  u) 

-00  -00 

OO  00 

f       Jh^2^y^  '^y     /   ^2^""^   '^"'  f7{t)   3(t  +  r      -  y  -  v) 

-OD  -00 

00  00 

4-     J  ivi(x)  dx    J  h^2^y)  dy  fpPTT^TPTrTT"^) 


-00  -00 

<©  00 


^       /h^l(x)   dx      I  h^2^y)  dy     j  h^iv)   dv 

-OO  -00  -OD 

(•)      f^Ct')    set'  +    f    +    X    -   y    -   v) 

00  00  00 

+     Jki^x)  dx    j  h^2^y^  ^    /  ^1^^^  d^ 

-<»  -00  -00 

(')   J^(t')   fgit'  +    t    +    X  -  y+  u) 


39 

CD 


-00  -00  -00  -00 

(•)  3(t')  5(1'+  f  +x-y+u-v)   (82) 

Now,  if  f4,  f5,  f5,  fy,  fs'  ^^  ^9'  ^^®  uncorrelated  noise 
functions,  and  satisfy  (22) 


CD 

Jf(t)     dt       = 


(83) 


then  the  cross  correlations  in  the  f d rst   eight  terms  of 
Equation    (82)  are  zero.     This   illustrates  one  of  the  advan- 
tages of  the  cross  correlation  function,    the  ability  to 
eliminate  the  uncorrelated  noise  from  the  final  answer. 
Ecfuat i on    (82)    is  now 


-<5o  -<D  -flo 


)   du 


oo 

(•)      j  hgtv)   dv    0g3(T  +    X   -  y  +  u   -  v)    (84) 

-00 

The  cross  power  spectrum  density  of  ^0102^'^^  "^"^  ^® 
obtained  by  using  Equation    (66).     I-^ltiplying  both   sides  of 
Equation    (84)  by  exp(-j(JT)  dT  f2TC  and  integrating  over 
T   gives 

<^DlD2<t^'    =     ^      f0^^^,iT).^l-iUT)ar     (85) 
-oo 

Using  the  variable 

^      =      T+x-y-j-u-v  (85) 
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in  Equation  (85)  gives 


CD 

-00 
00  00 

)       f  h^2(v)   exp(-jUy)  dv      J    ^^^(u)   expijUn)  du 


00 

( 

-do 


00  QD 

(.)       r  hgCv)   exp(.j^^v)   dv   2^     J    0ss(^'   exp(-j(JM)   dfX 
-00  -00  (87) 

Using  Equations    (66  and  68),   Equation    (87)  becomes 

^D1D2"^>     =     '4'^'  \2^^^     «I<^'  «2<'^'   ^ss»^Jg3, 

The  cross  power  spectrum  of  two  detectors  depends  on  the 
transfer  functions  of  the  two  measuring  systems,  the  transfer 
function  of  the  nuclear  system,  and  the  auto  power  spectmm 
of  the  source. 

Moore  (38.)  has  independently  developed  a  general  equa- 
tion for  the  auto  power  spectra  of  a  noise  field  that  is 
similar  to  Equation  (88).  T^oore,  however,  ignores  the  ef- 
fects of  the  measuring  system. 

Repeating  the  same  development  as  that  above  for 
^D1D3^^^  gives 

^D1D3«^'  =  «:;i"^'  \3'^'  '^I'e^'  «3'^'  ^as''^'  ^''^ 

If  Equation  (89)  is  now  divided  by  Equation  (88),  the  result 
is 
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Odid3^^^  _  HJ^^CJ)  \,iU)  H^((J)  n^iU)  j)33((J) 

^D1D2^^^     <l(^^  \2^^^  ^1^^^  ^2^^>  Cs^^^ 

(90) 

It  is  clear,  if  the  third  measuring  system  has  the  same 
transfer  function  as  the  second,  and  if  ^ss^^^  ^^  ^^® 
same  for  both  cross  power  spectrum  measurements,  that 
Equation  (90)  reduces  to 


^0102^^5     H2(U) 


=    HggCCJ)  (91) 


^ere 


H32((J) 


G(r2,r^,U) 
Aaj.(r2,r^,(J) 


(92) 


the  same  as  Equation  (72). 

Equation  (91)  indicates  that  the  transfer  function  of 
an  unknown  system  may  be  accurately  determined  by  the  ratio 
of  the  cross  porter  spectra  between  detectors. 

If  0D1D2^'^^  ^^^  0D1D3^'^^  ^^^  ^®  measured  at  the 
same  time,  there  is  no  restriction  on  the  source  behavior 
or  on  the  measuring  system.  In  actuail  practice,  it  is 
highly  desirable  to  have  the  band  width  of  the  measuring 
system  wider  than  the  band  width  of  the  unknown  system,  but 
this  is  not  critical  when  using  the  above  approach. 
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If  the  two  cross  correlations  cannot  be  performed  at 
the  same  time,  it  is  necessary  to  use  a  pseudo-random  source 
to  insure  that  it  repeats  the  same  input  for  both  measure- 
ments. The  frequency  content  of  the  source  should  be  wide 
enough  to  insure  that  the  entire  band  width  of  the  unknown 
system  is  excited. 

It  is  interesting  to  note  that  if  the  band  width  of 
the  measuring  system  is  wider  than  the  band  width  of  the  un- 
known system,  and  if  the  source  is  a  white  noise  over  the 
band  width  of  the  unknown  system,  both  ^gg(CJ)  and  ^j^^U) 
are  effectively  constant.  Then  Equation  (88)  becomes 

$D1D2^^^  =  ^ml  %^^^  ^2^^^  ^ss 

=  k'|h^ (1^)12  Hg^CU)  (93) 

Any  phase  shift  in  ^^ij)2^U)   is  caused  only  by  H2i(CJ), 
while  the  attenuation  of  <|>did2^^^  depends  on  E^^iLD^^^U) 
Even  if  the  transfer  function  of  the  measuring  system  is  not 
flat,  nor  the  source  white,  the  phase  shift  in  $diD2^^^  ^^ 
still  caused  only  by  U^i^iU), 

Since  the  above  assumptions  about  the  measuring  system 
and  source  may  not  be  valid,  it  is  usually  better  to  use 
Equation  (91),  even  with  the  requirement  of  an  additional 
measurement . 

Auto  power  spectrum  of  one  output 

The  auto  power  spectrum  of  Detector  1  is  of  interest. 
If  4)0102^^^  ^"^  $D1D3^^^  ^^®  measured  at  two  different 
times,  (|)r)in-i  ((^)  niay  be  checked  to  insure  that  the  fre- 
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quency  content  of  the  source  was  the  same  for  both  measure- 
ments. If  $gg(^)  is  different,  $didi(^)  may  be  used  as 
a  reference  in  comparing  (J)  ^^^^2^^  J  ^^^  $D1D3^^^' 

The  relation  between  <Pq]^qi(Cu' )  and  ^ggC^)  may  be 
obtained  using  the  same  procedure  as  used  in  the  previous 
section.  The  definition  (from  Equation  (64))  for  the  auto 
correlation  is 

T 

^DlDl(T)  =  T^X  2^  /^lO^tJ  flO^t  +  T)  dt    (94) 

-T 

Substituting  Equation    (76)   for  f2o(t),   or  better  yet,   using 
Equation    (82)   and  replacing  all  the  5,    8,    and  2  subscripts 
by  4,    7,   and  1,    respectively,   gives  0D1D1^^^»     ^^  ^^^  ^^® 
cross  correlation  functions  between  the  various  noise  func- 
tions and  the  source  are  equal  to  asero    (1^),   then 


^nim(r)    =   0.^(r)  + 


DlDl^  I    /  T-yy 


(•)  044(r  +  X  -  y: 


00  -OD 

OD  00 


L-OD 


GO  OD 

(•)       r  h3^(u)  du       r  ^i(v)  dv  0gg(r  +     X  -  y  4-  u   - 


v) 


(95) 

To   obtain  the  auto  power  spectrum,   both  sides  of  Equation 
(95)   are  multiplied  by  exp{-jU  T )  dT/27T  and  integrated 
over     I     (Equation    (66)).     Then,   using  the   same  reasoning 
presented  in  the  previous   section    (Equations  (85,    86,    87, 
and  88)), 
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ft  n    =  ^     ((. 
DlDl 


If  the  noise  inputs  are  white,  their  power  spectra  are  con- 
stant (22.).  Thus 

+  |«ml<^'|^  |Hi(U)|2  (J)33(U)   (97) 

So,  even  if  cj^gs^^^  cannot  be  checked  directly,  $didi{^) 
will  be  the  same  for  each  ran  if  (pggdJ)  is  also  the  same. 


CHAPTER  III 
EXPERIMENTAL  METHOD  AND  ANALYSIS 

The  analysis,  procedures,  and  considerations  nec- 
essary for  the  experimental  verification  of  the  spatially 
dependent  transfer  function  equations  developed  in  Chap- 
ter II  are  presented  and  discussed  in  this  chapter.  'Ihe 
method  used  to  obtain  the  experimental  data  is  discussed 
first,  followed  by  a  description  of  the  data  analysis, 
The   experimental  equipment  is  described  in  Chapter  IV. 

Experimental  Method 

The  experimental  data  was  obtained  using  the  method 
outlined  in  the  statistical  theory  section,  in  Figure  1 
and  in  the  development  of  Equations  (88  and  90). 

Ihe  flow  of  information  through  the  experimental 
system  is  shown  in  Figure  2.  Ihe  shift  register  gen- 
erates a  pseudo -random  binary  signal,  which  is  used  to 
turn  the  neutron  generator  on  and  off,  producing  the  in- 
put signal.  2ie  input  signal  travels  through  the  nuclear 
system  to  the  neutron  detectors.  The  neutron  detection 
system  is  a  pulse  system  vihich  produces  an  output  pulse 
of  constant  height  for  each  neutron  detected.  The   pulse 
rate  is  measured  by  a  fast  response  count  rate  circuit 
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v^ich  prodaces  an  analog  voltage  proportional  to  the  neu- 
tron density  at  the  detector.  This  signal  is  conditioned 
by  a  wide-band  band  pass  filter  to  remove  the  dc  compon- 
ent and  the  components  higher  in  frequency  than  one  half 
the  sanpling  rate  of  the  data  acquisition  system.  The 
conditioned  signal  is  sampled  by  the  data  acquisition  sys- 
tem, and  the  san^^les  are  stored  on  magnetic  tape  in  digi- 
tal form.  The  data  is  processed  by  an  IBM  709  computer 
at  the  University  of  Florida  Computing  Center,  using  the 
procedure  discussed  in  the  second  section  of  this  chapter. 

Jnpv^t  giqnal 

Equations  (89  and  90)  show  that  almost  any  input 
signal  may  be  used  to  measure  the  experimental  transfer 
function,  as  long  as  it  has  the  desired  band  width.  The 
input  actually  used  is  discussed  below. 

'Jhe  input  to  the  nuclear  system  is  controlled  by  a 
predetermined  pseudo- random  binary  signal  which  turns  the 
neutron  source  on  and  off.  'iixe  signal  is  generated  by  an 
eight  stage  shift  register,  with  the  correct  feedback  cir- 
cuitry for  producing  a  maximum  length  sequence  of  255  At 
intervals. 

'She   maximum  length  sequence  generated  by  a  shift 
register  has  a  special  auto  correlation  function,  which 
is  shown  in  Figure  3.  '-Qie  integration  required  to  cal- 
culate the  auto  correlation  function  (Equation  (54))  is 
performed  over  an  integral  number  of  periods,  T^,  of  the 
shift  register.  'Ihe  function  0gg  (  T )  has  a  constant 
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Figure  3.  Typical  maximum  length  sequence  and  its  auto 
correlation  function 
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negative  conponent  of  -1/N  magnitude,   N  being  the  number 
of   intervals   in  the  sequence.     With  every  Tj^,    the  period 
of  the  sequence,    the  shape  of  0gg(T)   approaches  a  delta 
function.      Ihis   characteristic  of  0gg(T')  has  been  used 
iX§J   ^o  measure  the  impulse  response  of  a  reactor  by  cross 
correlation  of  the  input  and  the  output.      If  0gg(T   -    V) 
in  Equation    (70)   is   replaced  by  a  delta  function,    the  cros£ 
correlation  becomes 

00 

0io(r)  =|h(t)  6(T  - 1)  dt 

-00 

=  h(r)At  (98; 

if  the  settling  time  of  h(T)  is  less  than  T^^.  The  im- 
pulse response  is  obtained  as  a  function  of  the  delay 
time. 

Iliis  particular  property  of  the  input  signal  is  not 
used  directly  in  this  study  for  several  reasons.  First, 
the  determination  of  the  spatially  dependent  transfer 
function  requires  accurate  measurements  to  show  the  phase 
shift  and  attenuation  of  the  neutron  disturbance  as  it 
passes  the  detectors.  Thus,  the  approximation  of  0gg(T) 
by  a  delta  function  is  not  an  accurate  one  for  this  vrork. 

Second,  the  system  being  measured,  represented  by 
h(T),  must  have  a  band  width  narrower  than  0gg(T)  if 
Equation  (98)  is  to  be  a  good  approximation  of  Equation 
(70).  Another  way  of  stating  this  is  that  the  time  con- 
stant of  the  system  must  be  less  than  the  period  of  the 
input  sequence,  usually  at  least  five  times  less. 
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Third,   h(T)    should  not   change  significantly  over  one 
At   int-erval,   the  average  width  of  the  triangle  on  the 
©ss^*^^   curve.      This  means  that  the  time  constant   of  the 
system  must  be  several  times  larger  than  At, 

The  auto  power  spectrum  of  the  source,    0ss^CJ) 
(Figure  4),   has   several   interesting  features.      If   0      (T) 
is  treated  as  an  aperiodic  function,    viiich   is  valid  if  the 
settling  time  of  the  system  is  less  than  the  period  of  the 
input   sequence  T^^,    the  spectrum  is  continuous  and  is    (using 
Equation    (66)) 


27r 


sin 


(J  At 


uA 


(99) 


However,  if  the  settling  time  of  the  system  is  greater 
than  T^,    0gg(T)  is  periodic,  and  the  auto  power  spectrum 
is  discrete  (using  the  periodic  Fourier  expansion.  Equa- 
tion (A2)) 


^ss^nU^^ 


2 


2 


(100) 


The  fundamental  frequency,  CJt,  corresponds  to  T,.  The 
locus  of  the  discrete  spectrum  components  has  the  same 
shape  as  Equation  (99),  the  continuous  case.  The  half 
power  frequency  is  approximately  0,44 /At  for  both  cases. 
It  is  clear  that  the  period  of  the  input  sequence  T^  should 
be  larger  than  the  settling  time  of  the  system  if  at  all 
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possible,  so  that  the  input  spectrum  will  appear  continuous 
to  the  system  being  investigated. 

Sampling  theory 

The  perturbations  in  the  neutron  density  caused  by 
the  input  signal  and  uncorrelated  noise  are  presented  in 
analog  form  by  the  measuring  system  and  sampled  at  a  high 
rate  by  the  data  acquisition  system. 

Since  the  analog  signal  is  sampled  a  finite  number  of 
times,  it  cannot  be  completely  reconstructed.  In  terms  of 
frequency,  the  reconstructed  signal  cannot  have  a  frequency 
content  hii^er  than  one  half  the  sampling  rate.  On  the 
other  hand,  if  the  analog  signal  has  frequencies  greater 
than  one  half  the  sampling  rate,  they  \-nll   contribute  false 
information  on  the  amplitudes  of  lower  frequencies  through 
the  process  known  as  aliasing  (28,).  Thus,  the  sampling  rate 
should  be  at  least  twice  as  high  as  the  band -pass  of  the  un- 
known system  to  obtain  the  maximum  amount  of  information. 
If  this  condition  ns  not  met,  because  of  limitations  in  the 
sanqDling  rate  or  lack  of  interest  concerning  the  high  fre- 
quency response,  the  input  to  the  sampling  equipment  must  be 
filtered  to  eliminate  the  frequency  components  above  one 
half  the  sampling  rate. 

If  it  is  desired  to  process  data  already  recorded,  but 
only  the  low  frequency  components  are  of  interest,  it  is 
possible  to  skip  samples,  using  only  every  n^"  sanple.   Now, 
however,  it  might  be  necessary  to  numerically  filter  the 
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data  to  remove  any  high  frequency  components  present  above 
one  half  the  new,  or  effective,  sampling  rate. 

The  variance  of  the  sanqsle  mean  can  be  used  to  esti- 
mate the  minimim  number  of  samples  needed  to  have  the  cor- 
rect function  average.  Lee  iZZ)   shovs 


^2 

^1 


=  ^Cl  (101) 


where  Cfl  is  the  variance  of  the  sample  mean,  (f  ^  is   the 
variance  of  the  function  being  sampled,  and  n  is  the  number 
of  samples.  Thus,  as  a  large  number  of  samples  is  taken. 


0'J  approaches  zero.  The  error  in  ^  may  be  taken  as Vl/^- 
o 

Data  Analysis 

The  data  stored  on  magnetic  tape  by  the  data  acquisi- 
tion system  corresponds  to  f^Ct)  and  f2("t),  the  two  output 
functions,  in  Equation  (64).  The  equations  used  to  numeri- 
cally process  this  data  and  compare  the  experimental  and 
theoretical  spatially  dependent  transfer  functions  are  dis- 
cussed in  this  section. 

Experimental  correlation  functions 

The  equations  developed  in  the  statistical  theory  sec- 
tion cannot  be  used  directly  in  a  practical  experiment,  since 
the  data  is  processed  digitally,  and  the  integration  times 
must  be  finite. 

Integration  time  errors.  Equation  (79)  may  be  approx- 
imated by 
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0D1D2<'7')   -  r/^lO^tJ  fii(t  +r)  dt 


(102) 


for  a  finite  integration  time  T. 

The  integration  time  necessary  to  obtain  a  reasonably 
accurate  answer  for  <^did2^'^^  ^^^  ^^^"^   investigated  by  oth. 
ers.  Balcomb  (15)  treated  the  cross  correlator  as  a  band 
pass  filter,  and  showed  that  the  "improvement  factor"  is  at 
least 


r-q/Ni     V  L 
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IF  is  defined  as  the  ratio  of  the  signal  to  noise  ratio  at 
the  input  to  the  same  ratio  at  the  output  of  the  cross  cor- 
relator, T  is  the  integration  time,  and  L  is  the  settling 
time  of  the  system.  The  above  derivation  is  based  on  a 
cross  correlator  designed  to  utilize  Equation  (98).  By 
making  T  =»  L,  the  uncorrelated  noise  is  averaged  over  a 
longer  period  than  the  correlated  signal.  Since  the  un- 
correlated noise  is  assumed  to  satisfy  Equation  (83),  its 
contribution  to  the  answer  is  decreased  with  greater  T. 
^ifhile   Equation  (103)  defines  an  improvement  factor,  it 
gives  no  information  on  the  amount  of  noise  present  before 
the  signal  enters  the  cross  correlator. 

Since  the  input  signal  discussed  earlier  is  a  periodic 
pseudo-random  function,  no  additional  information  is  gained 
about  the  system  by  integrating  for  times  longer  than  one 
period  of  the  input.  However,  the  3/N  ratio  may  be  sub- 
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stantially  improved  by  averaging  the  uncorrelated  noise  over 
several  periods, 

T?ajagopal  ilTj   approached  the  error  problem  by  con- 
sidering the  ratio  of  the  cross  correlation  function  for  an 
infinite  T  to  the  square  root  of  its  variance  for  a  finite 
T.  Considering  a  simple  lag  filter  as  the  unknoT-m  system 
T^rlth  a  white  noise  input,  the  signal  to  noise  ratio  from  the 
cross  correlator,  as  defined  in  Reference  (17.),  is 


crf^.^cr.TJ     ^|[i  +  e.p(-«r)] 


for  T=»=»T.  Thus,  as  in  Equation  (103 J,  the  3/N  is  pro- 
portional to  the  square  root  of  the  integration  time  T. 

As  the  delay  time  f  approaches  the  time  constant  of  the 
lag  system,  1/a,  the  integration  time  T  becomes  very  large 
for  what  seems  to  be  a  reasonable  error  limit.  For  example, 
for  S/N  equal  to  100,  T  equals  10^  seconds  for  af  =  0.625, 
For  S/N  equal  to  2,  T  becomes  a  reasonable  43  seconds.  Thus, 
even  though  the  cross  correlation  function  is  still  quite 
large  (exp (-0.625)  =  0.535),  the  error,  as  defined  in  Equa- 
tion (104),  con  also  be  quite  large  for  reasonable  integra- 
tion times.  Equation  (104)  is  specialissed,  however.  For 
inputs  that  have  a  finite  band  width,  the  error  will  be 
much  less  than  that  calculated  by  Equation  (104). 

Stem  (^)  presents  a  detailed  analysis  of  three 
sources  of  error  for  the  cross  correlation  function  as  re- 
presented by  Equation  (70).  This  analysis  is  applicable  to 
the  cross  correlation  calculations  performed  herein. 
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Defining  the  relative  error  as  (^) 

0io'  -^o' 


R     = 


(105) 


Stem    (21)   shows  that,    for  Ti/N  «=  <  T^-c -c  Tt^, 


r2  =  :^ 

2T 
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fT        2TT  7/2      [ 
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o-^R 


(106) 


To  is  the  time  constant  of  the  unknown  system,   and  T^  is 
the  period  of  the  pseudo-random  input   signal.     The  number 
of  At  intervals  in  the  input   signal  is  N,    1]  is  the  duty 
cycle  of  the  input    (-0.5),   and  t  is  the  average  counting 
rate.     The  average  source  strength  is  3^,    and 


A    = 


v^  -  V 


4-1 


(107) 


where  V   is  the  number  of  neutrons  emitted  per  fission  and 
p  is  the  reactivity  of  the  nuclear  system. 

The  terms  on  the  right  in  Equation  (106)  are  due  to 
the  finite  integration  time  T,  the  finite  counting  rate, 
and  the  natural  neutron  population  fluctuation,  or  the  un- 
correlated  noise  (21).  Ecfuation  (106)  may  be  used  to  esti. 
mate  the  required  integration  time  T  for  any  value  of  R. 

As  an  example  a  medium  of  light  water  with  the  fol- 
lowing typical  values  may  be  considered. 

Tr  =   2  X  lO-'^sec    ri  =  0.5    ft     =      0 

r   =   5  X  10^  counts/sec    S^  =  5  x  10*7  neu/sec 


0.1  sec 


So  = 

255 
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Then 

9     10-'^   5.2  X  10-5   2.04  X  10-^ 
R^  =  +  ^ 

T        T  T 

=  3.56  X  lO-'^/T  (108) 

Therefore,  for  an  error  of  one  per  cent,  the  integration 
time  T  is  3.56  seconds. 

This  is  in  sharp  contrast  to  the  10  seconds  obtained 
from  Equation  (104)  for  T/Tr  equal  to  0.625.  However, 
there  are  several  differences  between  the  two  approaches, 
preventing  them  from  being  compared  directly.  Rajagopal 
(17)  used  an  input  having  an  infinitely  wide  spectrum, 
while  3tem  (^)  used  an  input  having  a  finite  band  width, 
Rajagopal  considered  only  the  finite  integration  time 
error,  for  variable  T,  vith  a  simple  lag  system.  Stem 
considered  three  error  sources,  but  for  a  fixed  Tq» 

Probably  the  greatest  difference  is  that  Stern  con- 
sidered a  periodic  pseudo-random  source.  As  shown  in 
Reference  (25.),  the  variance  of  the  correlator  output  is 
smaller  if  a  periodic  pseudo-random  binary  input  signal  is 
used,  rather  than  a  purely  random  input.  The  variance  is 
reduced  even  more  if  the  integration  time  is  equal  to  an 
integral  number  of  periods  of  the  input.  Thus  it  would 
appear  that  the  purely  random  input  is  not  the  best  input 
signal  for  experimental  work.  The  infinitely  wide  band 
width  adds  unnecessary  noise,  requiring  greater  integration 
times  to  reduce  the  error. 
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It  is  clear  that  the  error  analysis  of  cross  cor- 
relation functions  is  not  an  exact  procedure.  However, 
consideration  of  the  discussion  presented  above,  parti- 
cularly that  of  Stem  (23.),  enables  a  reasonable  estimate 
of  the  error  to  be  made. 

Digital  analysis.  The  actual  calculation  of  the  cross 
correlation  function  was  done  numerically,  using  the  Uni- 
versity of  Florida  IBM  709  computer.  The  sampled  output 
functions  from  the  unknovm  system  were  read  into  the  com- 
puter and  the  cross  correlation  function  calculated  at 
discrete  delay  times  by  (from  Ecpiation  (102))  (31) 

I 
0DlD2^^Ar)  =  A^^fiodAt)  fiidAt  +  nAT) 
^~^  n  =  0  to  N 

I 
=  A^f^^(iAt)  f;^Q(iAt  +  nAT) 

^~  n  =   0  to  -N  (109) 

The  integration  time  T  is  l^t,   the  delay  time  is  nA  T  , 
and  A  is  a  normalizing  factor.  The  delay  time  interval, 
AT#  is  equal  to  A^  ^^^  is  determined  by  the  sampling 
rate  of  the  data  acquisition  system.  N  is  usually  much 
less  than  I.  The  functions  0did2^T)  ^^'^^diqs^T  )   are 
calculated  the  same  way,  using  the  appropriate  output  func- 
tions in  Equation  (109), 

Experimental  Power  Spectra 

The  power  spectra  are  calculated  by  performing  a 
numerical  Fourier  transformation  of  the  correlation 
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functions  (from  Equation  (66)).  Balcomb  (16.)  investigated 
the  effects  of  several  interpolation  schemes  on  the  Laplace 
transform  of  a  function  having  the  form  f  (n^x)  =  exp(-an^x) 
The  result  was  compared  with  the  known  transform  of  exp(-ax). 
His  conclusion  was  that  a  linear  interpolation  was  as  ac- 
curate as  higher  order  polynomial  interpolations,  and  was 
much  superior  to  no  interpolation  of  f  (nAx)  at  all. 

Since  0Q2D2^'^  ^  ^^^   expected  to  have  a  functional 
form 


fit)    =      ^A.  exp(-ajr|)  (110) 

i=l 

it  was  decided  to  use  a  linear  interpolation  of  (piT). 
Therefore,  the  numerical  form  of  the  Fourier  transform 
becomes 


-Tn)] 


exp(-j(jr)   dT      (111) 


where  0n=  ^D^^'^n^  and  AT  =  T^^^  -  T^  =  constant. 
Performing  the  integration  and  simplifying  the  sum  results 
in 


<|oiD2(U)     2     exp(.JUr.N) 


•J0.II       0-II+1   -0 

u      '        AT 
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N-1 
+   ;    exp( 
n=-N+l 


+  exp(-j 


(112) 


^j)203(U)  and  ^  DiT)l  (U )  are  calculated  by  siibstituting 
the  appropriate  0  (T)  i'^  Equation  (112).  A  Fortran  sub- 
routine TRNS  was  programned  to  calculate  the  power  spectrum 
using  Equation  (112).  It  was  included  as  part  of  the  gen- 
eral program  (31)  for  processing  the  experimental  data. 
See  Appendix  D  for  a  description  of  TRNS. 

Experimental  t.ransfer  .function 

The  transfer  function  is  calculated  using  a  modified 
form  of  Equation  (91).   Although  Equation  (91)  is  derived 
using  three  outputs  from  the  system  (Figure  1),  the  actual 
experimental  setup  (Figure  2)  contains  only  two  detection 
systems.   By  making  two  separate  recordings,  with  a  detector 
first  in  the  D2  position  and  then  in  the  D3  position, 
H23((ju/)  may  be  calculated  using  ^Pqilqi  (CJ)  as  a  normalizing 
factor. 


'izs"*^'  = 


r<l>DlD3«^'l 

L<i>DlDl((^'J 

Run   2 

Pd1D2«^'1 

L$D1D1<W)J 

Run  1 

(91a) 


The  numerical  errors  that  accumulate  as  the  data  is 
processed  cause  H2^{U)   to  become  less  accurate  as  the 
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power  spectra  are  attenuated  at  the  higher  frequencies. 
However,  as  noted  t^th  Equation  (93),  the  phase  shift  of 
the  system's  transfer  function  can  be  found  directly  from 
<I>D1D2(CJ)  or^DiD3(U). 

Ga,lcu3.ation  of  dynamic  parameters 

The  determination  of  the  dynamic  reactor  parameters 
was  aided  by  the  use  of  a  Fortran  code,  "General  Least 
Squares  Program",  originally  written  at  Oak  Ridge  National 
Laboratory  by  Bussing  and  extensively  modified  since  by 
Ccckrell  (J£)  for  use  by  the  Nuclear  Engineering  Depart- 
ment, University  of  Florida. 

Input  to  the  program  is  the  experimental  data  and  the 
theoretical  function  that  is  to  be  fitted.  The  program  is 
informed  which  parameters  in  the  theoretical  function  are 
to  be  varied  to  minimize  the  variance  between  the  data  and 
the  function. 

As  the  theoretical  function.  Equation  (57)  was  incor- 
porated into  a  main  subroutine  CALC,  and  several  minor 
subroutines,  C1?T,  GUEST,  GUESS,  TERM,  ERROR,  and  FINAL. 
The  computer  codes  are  discussed  in  detail  in  Appendix  D, 

It  is  not  possible  to  determine  the  five  parameters 
J?g,  J?^,  L^,  T,  and  k^^  simultaneously.  The  dependence 
of  the  transfer  function  on  them  is  such  that  no  single 
combination  of  4,  S-r,    L^,  T,  and  k^^  exists  for  mini- 
mizing the  variance.  However,  if  two  or  three  of  the 
parameters  are  known  from  other  sources,  the  remaining 
ones  may  be  estimated  using  the  General  Least  Squares 
program. 


CHAPTER  IV 
EXPERIMENTAL  EQUIPMENT 

The  equipment  used  to  generate  the  input  signal  and 
to  measure  the  output  of  the  subcritical  assembly  is 
described  in  this  chapter.  A  schematic  of  the  measuring 
system  is  shown  in  Figure  2. 

The  experiments  using  the  light  water  subcritical 
assembly  were  conducted  at  a  different  location  than  the 
experiments  using  the  heavy  water  subcritical  assembly. 
Thus,  some  items  of  equipment  were  different  for  the  t^vo 
experiments. 

Input  System 

The  input  system  contains  two  major  components:   the 
shift  register  and  the  neutron  generator. 

Shift  Register 

The  pseudo-random  binaiy  input  signal  is  generated 
by  a  shift  register  with  "modulo  2*  feedback.  By  using 
the  proper  feedback  logic,  the  maximum  length  sequence 
can  be  obtained.  This  sequence  has  the  special  auto 
correlation  function  described  in  Chapter  III.  Peterson 
(26)  discusses  the  correct  feedback  circuitry  for  obtain- 
ing the  maximum  length  sequence  for  shift  registers 
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containing  up  to  34  stages.  The  number  of  intervals,  N, 
contained  in  a  maximum  length  sequence  generated  by  a 
shift  register  having  m  stages  is 

N  =  2^-1  (113) 

References  (^3^  25,)  describe  in  some  detail  the  charac- 
teristics of  a  shift  register  used  to  generate  a  pseudo- 
random signal. 

The  shift  register  used  in  the  present  study  has  eight 
stages,  producing  a  sequence  having  255  intervals  before 
repeating.  A  schematic  of  the  shift  register  is  shown  in 
Figure  5.  The  /\t  time  interval  is  controlled  by  a  free 
running  multivibrator,  whose  speed  may  be  varied  from  1  cps 
to  10  KC. 

Neutron  Generator 

A  Norelco  neutron  generator  was  used  for  the  experi- 
ments with  light  water,  and  a  Texas  Nuclear  Corporation 
neutron  generator  was  used  for  the  heavy  water  experi- 
ment s . 

Norelco  neutron  generator.  The  Norelco  neutron  gen- 
erator (3^)  produces  neutrons  of  14  IffiV  energy  by  the  d-t 

8 
reaction,  with  an  average  source  strength  of  10  neutrons 

per  second. 

l!he  neutron  source  tube  contains  a  replenisher  sys- 
tem which  emits  hydrogen  isotopes  y,jb.en   it  is  heated  by  an 
electric  current.  The  isotopes  are  ionized  in  a  Penning 
ion  source  and  moved  to  the  one  stage  accelerating  system 
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by  the  ion  source  voltage.  The  accelerating  system  imparts 
enough  energy  to  the  ions  (125  KV)  to  overcome  the  electro- 
static potential  of  the  deuterium  and  tritium  nuclei  al- 
ready embedded  in  the  titanium  target,  thus  causing  the 
nuclear  reaction.  The  size  of  the  source  tube  is  3'"  x  20". 

The  neutron  generator  may  be  pulsed  by  controlling 
the  ion  source  voltage,  with  a  maximum  pulsing  rate  of 
about  3  KC  at  a  50  per  cent  duty  cycle.  The  signal  gen- 
erated by  the  shift  register  is  used  to  control  the  ion 
source  voltage  through  an  intermediate  voltage  matching 
amplifier. 

A  high  voltage  transformer  that  supplies  the  accel- 
erating voltage  is  located  close  to  the  source  tube.   A 
control  rack,  containing  the  electronic  equipment  and  ion 
source  voltage  transformers,  is  connected  to  the  source 
tube  by  a  30-foot  cable.  This  length  allows  the  rack 
and  operator  to  be  located  a  safe  distance  from  the  neutron 
generator. 

Texas  Nuclear  Corporation  neutron  generator.  The  TNC 
neutron  generator  produces  neutrons  by  either  the  d-t  or 
d-d  reactions.  The  d-d  reaction  was  used  for  the  heavy 
water  experiments. 

While  the  TNC  neutron  generator  is  physically  differ- 
ent from  the  Norelco  neutron  generator,  the  sequence  of 
events  required  to  produce  neutrons  is  the  same. 

The  replenisher  system  is  a  pressurized  tank  con- 
taining deuterium  gas.  The  deuterium  is  ionizing  in  a 
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Penning  ion  source  and  is  moved  to  the  acceleration  system 
by  the  ion  source  voltage.  The  high  energy  ions  leave  the 
acceleration  system  and  strike  the  deuterium,  titanium 
target,  producing  neutrons. 

The  shift  register  controls  the  pulsing  of  the  neu- 
tron generator  through  an  intermediate  circuit.  The  random 
width  pulse  from  the  shift  register  is  differentiated,  and 
on  and  off  pulses  are  supplied  to  the  neutron  generator 
to  control  the  beam  deflection  plates.  The  plates  are 
located  between  the  ion  source  and  the  acceleration  system. 
A  higher  pulsing  rate  is  possible  with  this  generator  than 
with  the  Norelco  generator. 

Nuclear  Systems 

A  light  water  subcritical  assembly  and  a  heavy  water 
subcritical  asseinbly,  each  containing  various  amounts  of 
natural  uranium,  were  used  as  the  test  media. 

Light  Water  Subcritical  Assembly 

The  light  water  subcritical  assembly  consists  of  a 
right  cylinder  tank  containing  water,  an  aluminum  support 
grid,  and  natural  uranium  fuel  rods.  The  main  features 
of  the  tank  are  shown  in  Figure  5. 

The  neutron  generator  was  installed  inside  a  3  1/4 
inch  horizontal  pipe  in  the  bottom  of  the  tank,  such  that 
the  neutron  source  was  on  the  longitudinal  axis  of  the 
assembly.  With  no  fuel  elements,  the  effective  radius  of 
the  tank  is  52  cm.  When  the  uranium  is  added,  the  core  is 
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Figure  6.      Light  water  subcritical  asseinbly 
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surrounded  by  a  water  reflector.  For  a  full  loading  of 
fuel,  the  effective  radius  of  the  core  is  37  cm.  The  core 
is  immediately  above  the  neutron  generator  and  extends  to 
the  surface  of  the  water,  142.9  cm  from  the  source. 

The  support  grid  holds  the  uranium  in  a  triangular 
lattice  4  cm  on  a  side.  Two  hundred  four  fuel  elements 
make  up  a  full  load,  with  one  diagonal  row  of  elements 
missing  to  leave  room  for  the  neutron  detectors.  The 
natural  uranium  slugs  are  housed  inside  1  inch  aluminum 
tubing,  which  is  held  upright  by  the  grid. 

The  volume  fractions  for  the  full  load  are:   %20  ~ 
0.42133,  F^L=  0.13058,  and  F^  =  0.41809. 

Heavy  Water  3ubcr;it-i  cal  Assembly 

The  heavy  water  subcritical  assembly  (Figure  7)  is 
essentially  the  same  as  the  light  water  subcritical  as- 
sembly. The  main  differences  are:   a.  The  target  of  the 
neutron  generator  is  placed  inside  a  2  1/2  inch  pipe 
through  the  center  of  the  tank.  b.  The  lattice  spacing 
is  12  cm,  resulting  in  volume  fractions  of:   pQgO"  0*9357, 
F;^L=  0.01784,  and  Fgp  =  0.04546,  for  a  full  loading  of  55 
fuel  elements,   c.  The  lattice  occupies  the  entire  volume 
of  the  DgO  tank.  Thus,  the  effective  radius  is  59  cm  for 
both  the  non -multi plying  and  multiplying  cases,  d.  The 
water  level  is  maintained  at  60.8  cm  above  the  source. 
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Figure  7.  Heavy  water  subcrjtical  assembly 
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Output  System 

The  variation  of  the  neutron  density  is  sensed  by  the 
output  system.  The  major  components  in  this  system  are  the 
detection  and  amplification  equipment,  the  count  rate  cir- 
cuits, the  band  pass  filters,  and  the  data  acquisition 
system. 

Neutron  Detectors 

The  neutron  detectors  were  placed  inside  aluminum 
conduit  to  protect  them  from  the  water.  The  conduit  was 
L-shaped,  supporting  the  detectors  in  a  horizontal  posi- 
tion above  the  neutron  generator.  The  height  of  the  in- 
dividual detectors  could  be  easily  changed  by  vertical 
adjustment  of  the  conduit. 

The  detectors  have  a  12  inch  active  length,  filled 
with  He^  to  4  atmospheres  of  pressure. 

Pulse  Amplifiers 

For  the  light  water  experiment,  preamplifiers  and 
linear  amplifiers  manufactured  by  Hamner  Electronics 
Inc.  were  used.  Preamplifiers  and  linear  amplifiers 
made  by  Nuclear  Chicago  were  used  for  the  heavy  water 
experiment.  Both  systems  supplied  constant  height,  one 
microsecond  long  voltage  pulses  to  the  count  rate  circuits 
at  a  rate  proportional  to  the  neutron  density  at  the  de- 
tectors. 
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Count  Rate  Circuit 

The  count  rate  circuit  was  specially  designed  for 
this  study.  The  specifications  were:   a,  a  fast  response 
circuit,  whose  output  would  duplicate  the  neutron  density 
variation  as  closely  as  possible;  b.  little  ripple,  or 
decay  of  the  output  voltage  between  input  pulses;  and 
c.  adequate  output  voltage  (2  to  10  volts). 

The  count  rate  circuit  has  five  counting  ranges, 
T/ijhich  are  shof«m  with  their  measured  break  frequencies  in 
Table  1.  Appendix  E  contains  detailed  information  on 
the  count  rate  circuit,  including  the  design  criteria, 
a  circuit  diagram,  and  a  plot  of  the  transfer  function  of 
the  circuit- 


TABLE  1 

Count  Rate  Meter 

Rangg. 

Cpuntsyse^. 

-3  4b  point 

1 

lOOK  to  250K 

5400  cps 

2 

40K:  to  lOOK 

2000  ops 

3 

20K  to  5 OK 

1000  cps 

4 

lOK  to  25K 

600  cps 

5 

4K  to  lOK 

350  cps 

FM  T^p^   ^e^ord^r, 

A  Minneapolis-Honeyr^rell  four-channel  tape  recorder 
was  used  for  storing  the  output  from  the  count  rate  cir- 
cuits during  the  heavy  water  experiments  and  for  trans- 
porting the  information  to  the  data  acquisition  system. 
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The  tape  recorder  has  a  band  width  of  0  to  10  KC 
at  a  recording  speed  of  30  inches  per  second;  a  speed  of 
1  7/8  inches  per  second  is  also  available.  An  internal 
calibration  circuit  and  compensation  channel  allows  the 
original  signal  to  be  reproduced  with  less  than  one  per 
cent  error. 

Band  Pass  Filters 

Krohn-Hite  ultra-low  frequency  band  pass  filters 
were  used  to  condition  the  analog  signal  from  the  count 
rate  circuit.  The  high  pass  filter  was  set  at  .02  cps 
to  remove  the  low  frequency  components  from  the  signal. 
The  low  pass  filter  was  set  to  remove  the  frequencies 
higher  than  one  half  the  sampling  rate  of  the  data 
acquisition  system. 

Data  Acquisition  System 

The  main  components  of  the  data  acquisition  system 
are  shown  in  Figure  8.  The  multiplexer  samples  two,  four, 
six,  or  eight  input  analog  signals  at  a  rate  of  15,000 
samples  per  second  in  the  low  density  tape  mode.  The 
twelve-bit  binary  analog  to  digital  converter  accepts  an 
input  from  each  channel  in  turn  of  -4.095  to  4.096  volts. 
The  format  unit  calculates  and  includes  in  the  final  bi- 
nary word  a  lateral  odd  parity  bit.  The  unit  also  cal- 
culates the  longitudinal  parity  check  sum  for  inclusion 
at  the  end  of  each  record,  and  performs  other  control 
functions. 
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The  resultant  computer  words  are  written  on  tape 
by  the  recorder  unit,  in  a  form  compatible  with  the  IBM 
709  computer.  The  system  can  also  produce  a  high  density 
tape  for  the  TB'l  7090  computer,  at  a  sampling  rate  of 
28,571  samples  per  second. 
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CHAPTER  V 

REsaurs 

The  results  from  the  theoretical  and  experimental 
study  of  the  spatially  dependent  transfer  function  of  a 
nuclear  system  are  presented  in   this  chapter. 

Theoretical  Results 

The  computer  code  SPAT,  described  in  Appendix  D, 
was  used  to  calculate  the  transfer  function  between  a 
thermal  neutron  detector  and  a  point  source  of  fast  neu- 
trons (Equation  (57))  for  four  cases,  presented  in  Figures 
9  through  13.  The  same  geometry  described  in  Chapter  II 
was  used,  with  values  of  c  =  150  cm,  and  R  =  50  cm.  The 
values  of  the  nuclear  parameters  used  in  the  calculations 
are  listed  in  Table  2. 

The  series  in  Equation  (57)  were  terminated  when  the 
magnitude  of  additional  terms  was  less  than  one  tenth  per 
cent  of  the  sum  of  the  absolute  values  of  the  preceeding 
terms.  The  curves  labeled  L  in  the  figures  represent  the 
lumped  parameter  transfer  function  for  the  nuclear  sys- 
tems. 
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TABLE  2 
Nuclear  Parameters  Used  for  the  Theoretical  Results 


Figure 

l\    cm2 

r,    cm2 

30 

f,.flsec 

200 

10 

^co 

9,   HgO^ 

8.12 

0 

10,   HgO^ 

8.12 

30 

200 

10 

0 

11,   DgO^ 

13,700 

110 

40,000 

46 

0 

12,   HgO  + 

Ur^ 

1.61 

40 

200 

10 

0.923 

13,   DgO  + 

Ur<= 

87.2 

133 

300 

40 

1.15 

b.  References  (li,  29.) 

c.  References  (30/  11/  li^ 

The  variation  of  the  phase  shift  with  frequency  and 
with  distance  from  the  source  is  of  interest.  At  100  cm 
from  the  source,  the  phase  shift  in  light  water  (Figure  9, 
Part  D,  and  Figure  10)  is  -720  degrees  at  800  cps  and  is 
-1440  degrees  at  2400  cps.  In  other  words,  four  complete 
cycles  of  the  2400  cps  neutron  wave  are  located  between 
the  source  and  the  detector.   As  the  detector  is  moved 
closer  to  the  source,  the  transfer  function  experiences 
a  maximum  negative  phase  shift  (Figure  9,  Part  C).  A 
point  is  reached  where  the  phase  shift  is  always  positive 
(Figure  9,  Part  B).  The  data  for  heavy  water.  Figure  11, 
also  shows  a  positive  phase  shift.  The  roultiplying  media 
data.  Figures  12  and  13,  shows  a  great  reduction  with  fre^ 
quenoy  in  the  amount  of  the  phase  shift. 
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Figure  9,  c  and  d.  Theoretical  transfer  function  for  light 
water 
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Figure  12,    a  and  b.      Theoretical  transfer  function  for  a 
light   water  subcritical  asseinbly 
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The  behavior  of  the  phase  shift  clearly  illustrates 
the  effect  of  frequency  on  the  wave  length  and  wave  velo- 
city. A  discussion  of  the  dependence  of  the  wave  velo- 
city and  wave  length  of  a  thermal  neutron  disturbance  on 
frequency  may  be  found  in  References  (5_,  6_) . 

A  comparison  of  the  spatially  dependent  transfer 
function,  3DTF,  and  the  lumped  parameter  transfer  function, 
LPTF,  shows,  in  general,  that  the  lumped  parameter  model 
is  not  an  adequate  representation  of  the  nuclear  system 
(see  Figures  9  through  13).  However,  it  appears  that  the 
two  models  can  give  similar  results  if  the  output  detector 
is  located  at  a  specific  distance  from  the  source. 

As  the  multiplication  factor  approaches  one,  the  two 
models  give  essentially  the  same  attenuation  and  phase 
shift  for  z  =  20  cm.  As  Y^    approaches  zero,  the  heavy 
water  case  can  still  be  represented  by  the  LPTF  for 
z  =  20  to  30  cm  and  for  frequencies  less  than  1000  cps, 
but  there  is  a  significant  difference  between  the  SDTF 
and  the  LPTF  for  the  light  water  case. 

Mclnerney  (40_)  has  shown  that  diffusion  theory  gives 
the  correct  spatial  dependence  for  fast  neutrons  when  the 
reciprocal  diffusion  length,  1/L,  is  approximately  equal 
to  the  total  cross  section,  2..^.  Mclnemey's  development 
helps  explain  why  the  LPTF,  which  is  derived  from  diffusion 
theory,  cannot  be  used  with  light  water  but  can  be  used 
in  a  light  water  moderated  reactor  at  specific  distances 
from  the  source.   As  uranium  is  added  to  the  water,  1/L 
approaches  L^,  fulfilling  the  requirement  stated  above. 


85 

In  using  Equation  (57)  to  determine  the  parameters 
of  a  nuclear  system  from  experimental  data,  it  is  of  in- 
terest to  know  the  number  of  terms  that  must  be  used  to 
adequately  fit  the  data.  An  indication  of  the  effects 
of  changes  in  z,  media,  and  frequency  on  the  number  of 
terms  was  obtained  by  considering  the  four  cases  listed 
in  Table  2.  The  results  presented  in  Table  3  were  ob- 
tained by  assuming  the  average  experimental  error  to  oe 
four  per  cent  and  using  this  percentage  as  the  accuracy 
criterion  in  3PAT. 

As  expected,  only  the  first  term  of  Equation  (57) 
is  needed  when  the  detector  is  at  a  sufficient  distance 
from  the  source.  The  higher  spatial  modes  are  attenuated 
more  than  the  fundamental  mode  and  can  be  neglected  at 
low  frequencies.  At  higher  frequencies,  though,  the 
hi^er  spatial  harmonics  are  not  attenuated  much  faster 
than  the  fundamental  mode.  They  must  be  included  in 
the  equation  even  for  large  z.  Figure  31  illustrates 
part  of  the  reason  for  this  effect. 

The  fundamental  root,  Q    ,   which  is  much  smaller  than 
the  other  roots  at  a  frequency  of  0  cps,  is  only  slightly 
smaller  at  2000  cps.  At  higher  frequencies,  at  least 
three  terms  in  the  i  series  must  be  included  in  the  equa- 
tion, and  possibly  several  terms  from  the  n  series  (see 
Equation  (57)).  The  spatial  modes  contributed  by  the  m 
series  are  usually  quickly  attenuated,  so  that  only  the 
fundamental  m  term  need  be  included. 
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TABLE  3 
Number  of  Terms  Required  for  Pour  Per  Cent  Error 
in  Equation  (57) 


Case 

Firegu^nqy, 

CPS 

0 

10 

100 

1000 

Figure  9,   H^O 

z     =     20  cm 

6 

7 

7 

7 

2     =     50  cm 

4 

4 

4 

5 

z     =  100  cm 

2 

3 

3 

3 

Figure  11,    D^O 

z     =     20   cm 

3 

3 

3 

5 

z     =     50  cm 

2 

2 

2 

3 

z     =  100   cm 

1 

1 

1 

2 

Figure  12,   H^O  + 

Ur 

z     =    20  cm 

3 

3 

3 

7 

z     =     50  cm 

2 

2 

2 

3 

2     =  100  cm 

1 

1 

1 

2 

Figure  13,    DgO  + 

Ur 

z     =    20  cm 

3 

3 

6 

9 

z     =     50  cm 

1 

1 

2 

2 

z     =  100  cm 

1 

1 

1 

2 
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As  the  multiplication  factor  approaches  one,  the 
fundamental  mode  becomes  dominant,  making  fewer  terms 
necessary  for  most  systems,  even  close  to  the  source. 

Unfortunately^  it  is  more  difficult  to  make  measure- 
ments in  that  portion  of  the  system  which  can  be  represented 
by  only  one  term  of  the  equation,  since  there  are  fewer 
neutrons  present  and  the  counting  rate  is  lower. 

Experimental  Results 

The  methods  used  to  check  the  accuracy  of  the  equip- 
ment and  the  computer  calculations  are  discussed  first. 
The  experimental  data  follows. 

Calibration 

Equation  (90)  shows  that  the  effects  of  the  measuring 
system  on  the  input  signals  are  cancelled  when  the  ratio 
of  the  cross  power  spectra  between  two  detectors  is  cal- 
culated. Thus  a  calibration  of  the  measuring  system  is 
not  necessary  in  the  usual  sense  of  the  word.  However, 
it  is  important  to  know  the  amount  of  noise  that  is  added 
to  the  input  signal  by  the  equipment  and  to  determine  the 
extent  of  the  numerical  errors  accumulated  in  the  data 
during  the  computer  processing. 

The  computer  code  was  checked  by  correlating  and 
transforming  known  mathematical  functions.  Once  the  pro- 
gram was  debugged,  the  error  between  the  calculated  and 
the  theoretical  results  (correlations  and  power  spectra) 
was  less  than  one  per  cent. 
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The  next  step  was  to  sample  the  output  from  the 
shift  register  with  the  data  acquisition  system,  and  pro- 
cess the  data  in  the  computer.  The  auto  correlation  of 
the  recorded  function  exhibited  the  characteristic  trian- 
gular spike  at  T  =  0  (see  Figure  3),  but  it  did  not 
remain  at  a  constant  -1/N  (-0.00392)  level  for  the  other 
delay  times.   Instead,  the  auto  correlation  function 
fluctuated  between  values  of  0.005  and  -0.02.  This  fluc- 
tuation is  caused  partly  by  numerical  error  in  the  calcu- 
lation and  partly  by  the  discrete  sampling  interval  of  the 
input  voltage.  The  majority  of  the  variation,  though, 
is  caused  by  noise  added  to  the  signal  by  the  electronic 
equipment.  The  equipment  involved  in  this  test  was  the 
shift  register,  an  intermediate  Schmidt  trigger,  and  the 
data  acquisition  system. 

Figure  14  shows  the  auto  power  spectrum  of  the  shift 
register  and  the  theoretical  curve  for  the  idealized  sig- 
nal. The  approximately  1  db  gain  in  the  20  to  200  cps 
frequency  range  is  caused  by  the  low  frequency  fluctuations 
noted  in  the  auto  correlation  function.  If  the  results 
presented  in  Figure  14  are  used  as  a  calibration  test,  it 
appears  that  -13  to  -15  db  is  the  approximate  limit  for 
accurate  calculation  of  the  power  spectrum. 

There  are  several  other  items  to  consider  when  judging 
the  accuracy  of  the  experimental  results.  The  number  of 
computational  steps  a  particular  sampled  function  has  gone 
through  affects  the  accuracy  of  the  function.  As  the 
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power  spectra  are  attenuated  at  the  higher  frequencies  and 
the  numbers  become  small,  the  error  is  increased.  Com- 
bining the  spectra  as  in  Equation  (91a)  can  produce  wildly 
fluctuating  results.  Since  the  functions  for  calculating 
(|)q^qj^((J)  and  4^DiD2^^^  ^^®  recorded  at  the  same  time 
and  then  subjected  to  the  same  computational  procedure, 
any  trends  introduced  by  the  process  tend  to  cancel  out 
when  <|)q3^q2^^^  ^^  divided  by  ^^^-^j^^iU)    in  Equation  (91a). 
This  process  is  particularly  helpful  when  the  correlation 
curves  are  discontinuous  at  tN  (see,  for  example.  Figure 
18),  as  the  discontinuity  causes  similar  bumps  on  both 
power  spectra  curves  (see  Figures  IS  and  20). 

Data 

The  correlation  functions  presented  in  this  section 
were  calculated  using  Equation  (109).  The  integration 
time,  lAt,  was  equivalent  to  ten  cycles  of  the  shift 
register.  Equation  (112)  was  used  to  calculate  the  power 
spectra.  The  time  span  of  the  integration,  -N  to  +N, 
was  varied  to  transform  only  that  portion  of  the  corre- 
lation data  which  was  of  interest. 

Runs  1  and  3  were  made  using  the  light  water  sub- 
critlcal  assembly  tank,  which  was  filled  with  water  only. 
The  cross  correlations  (Equations  (79,  80,  102,  109))  are 
shown  in  Figure  15.  \Vhen  only  an  "x"  appears,  it  coincides 
with  the  dot. 

The  correlations  0.  ( T  )  and  0to(T)  are  essentially 
the  same,  preventing  the  determination  of  the  spatial 
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effect  on  the  system's  transfer  function.  The  two  main 
reasons  for  their  similarity  are:  a.  the  rapid  decrease 
in  the  counting  rate  as  z  increased  prevented  a  sufficient 
Az  between  the  detectors,  and  b.  the  sampling  rate  of  the 
data  acquisition  system  is  not  rapid  enough  to  permit 
measurement  of  the  travel  time  of  the  neutron  disturbance 
between  the  detectors. 

Runs  2  and  6  were  made  in  the  heavy  water  subcritical 
assembly  tank,  which  contained  water  only.  Figure  16 
shows  the  cross  correlation  functions.  The  effect  of  moving 
Detector  2  to  Position  3  can  be  observed  in  the  0,4  milli- 
second shift  of  the  peak,  showing  that  measurement  of  the 
travel  time  of  the  neutron  disturbance  between  detectors 
is  possible. 

The  striking  feature  of  Figure  16  is  that  the  posi- 
tive correlation  peaks  appear  to  be  displaced  to  the  ri^t 
by  approximately  10  milliseconds.  The  reason  for  this 
displacement  is  discussed  later  with  the  heavy  water  sub- 
critical  data.  The  heavy  water  data  was  not  processed 
further. 

Figure  17  presents  the  data  taken  in  the  light  water 
subcritical  assembly  with  a  full  load  of  uranium.  The 
correlation  function  does  not  decay  to  zero  within  the 
limits  of  the  figure.  Even  for  much  longer  delay  times, 
03,2(T)  does  not  go  to  zero.  This  strange  behavior 
appears  to  be  caused  by  the  reflection  of  neutrons  from 
the  concrete  floor.  Other  structures  were  far  enough  away 
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from  the  tank  to  have  little  effect  on  the  neutron  dis- 
tribution in  the  tank. 

If  it  is  assumed  that  the  actual  source  variation, 
S(t),  is  the  intended  source  variation,  3^(t),  plus  a 
reflection  of  3^  from  the  floor,  the  actual  source  sup- 
plying neutrons  to  the  detectors  is 

S(t)  =  Si(t)  +  AS.  (t  -  T)  (114) 

The  intensity  of  the  reflected  source  relative  to  the  in- 
tended source  is  A,  and  T  is  the  time  lapse  before  the 
reflected  neutrons  re-enter  the  medium  and  pass  the  source 
position.  Using  the  same  procedure  as  that  used  in  Chap- 
ter II  to  develop  Equations  (88,  89),  the  cross  power 
spectrum  between  two  detectors  with  the  source  defined 
by  Equation  (114)  is 

(|),2((J)  =  <,((J)  ^^^iU)  nliU)  H^iU) 

(•)  (J)gg(U)[l  +  a2  +  2A  cos((jT)]    (115) 

The  bracketed  term  shows  the  effects  of  the  reflection. 
If  A  is  zero.  Equation  (115)  is  the  same  as  Equation  (88). 
The  cos((jT)  term  causes  a  'Treating"  of  the  power  spectra. 
The  transform  of  the  correlation  function  in  Figure  17 
produces  such  an  effect.  Resonance  peaks  occur  approxi- 
mately every  118  cps,  which  corresponds  to  the  8.4  milli- 
second delay  time  to  the  first  reflection  peak  in  Figure 
17.  Figure  17  actually  contains  other  reflection  peaks 
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at  multiples  of  8.4  milliseconds,  for  both  positive  and 
negative  delay  times.  A  single  value  for  A  could  not  be 
found  from  the  experimental  data.   A  appears  to  be  a 
function  of  frequency,  varying  from  0.2  at  200  cps  to 
0.4  at  450  cps. 

A  further  check  of  the  assumption  that  the  displaced 
correlation  peaks  result  from  reflections  can  be  made  by 
estimating  the  average  energy  of  the  reflected  neutrons. 

Using  the  definition  of  the  group  velocity  (41) 


^^  (116) 


""^     b[Real  (Bp)'J 
and  the  solution  for  the  buckling  from  simple  monoenergetic 


diffusion  theory. 
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The  average  neutron  energy  is  approximately 
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the  diffusion  coefficient  is  D,  3^  is  the  transverse 
buckling,  Zp  is  the  removal  cross  section,  and  m  is  the 
neutron  mass.  Since  the  reflected  neutrons  must  have  re- 
latively high  energies  to  survive  the  trip  into  the  floor 
and  back  to  the  detectors  through  the  water,  the  values  of 
D  and  E^  used  in  the  calculation  are  representative  of 
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fast  neutrons  and  are  0.2  cm  and  0.1  1/cm  (44).  The  group 
velocity  was  estimated  to  be  1.7  x  10^  cm/sec  by  considering 
the  time  delay  between  the  correlation  peaks  shown  in 
Figure  18  and  the  distance  between  Detectors  2  and  3. 
Using  these  values  in  Equation  (118)  gives  an  average  ener- 
gy of  the  neutrons  between  the  two  detectors  as  approxi- 
mately 0.19  ev.  IVhile  this  energy  is  low,  it  is  several 
times  greater  than  thermal  energy  and  it  does  indicate 
that  the  displaced  correlations  could  be  due  to  reflected 
neutrons.  Clearly,  further  study  of  the  reflections  of 
neutron  waves  and  groups  from  both  the  experimental  and 
theoretical  standpoint  is  needed.  Clarification  of  the 
mechanism  involved  in  the  spreading  and  reflecting  of 
neutron  groups  should  prove  most  valuable. 

Since  the  development  of  the  spatially  dependent 
transfer  function  included  the  assumption  that  no  neutrons 
re-entered  the  system,  the  cross  power  spectra  were  ob- 
tained by  transforming  that  part  of  the  correlation  curves 
that  is  close  to  T  =  0.  Unfortunately,  since  the  cor- 
relation curves  are  rarely  equal  to  zero  at  the  ±N  points, 
a  discontinuity  is  present  which  adds  a  small  beating 
effect  to  the  power  spectra. 

Figures  18  through  23  contain  the  data  obtained  from 
the  light  water  and  natural  uranium  subcritical  assembly 
used  to  determine  the  experimental  spatially  dependent 
transfer  function.  Only  that  portion  of  the  correlation 
curves  between  -4  to  +4  milliseconds  was  used  to  obtain 
the  power  spectra  curves.  In  Figures  19  and  20,  many  of 
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the  dots  coincide  v/ith  the  x's.  The  points  on  the  two 
figures  begin  to  become  scattered  for  frequencies  greater 
than  800  cps  and  for  an  attenuation  greater  than  15  db. 
Also,  the  discontinuities  in  the  cross  correlation  func- 
tion cause  bumps  on  the  spectra  curves,  starting  at  300 
cps.  The  correct  curve  is  an  average  of  these  bumps. 
The  data  beyond  300  cps  may  be  regarded  with  suspicion 
and  that  beyond  800  cps  may  be  disregarded. 

The  cross  power  spectrum  curve  for  Run  39,  where 
the  detectors  are  close  together,  has  essentially  the 
same  shape  as  the  auto  power  spectra  of  Detector  1  (see 
Figures  19  and  20).  Referring  to  Equations  (88,  97), 
the  above  fact  indicates  that  there  is  very  little  un- 
correlated  noise  in  the  system;  in  other  words,  Kyy  and 
K44  are  very  small.  Thus  (poiDi^^^  ^^V  ^®  safely  used 
in  Equation  (91a)  as  a  substitute  for  (DggdJ) . 

Figures  22  and  23  contain  the  results  from  Equation 
(91a).  The  attenuation  varies  smoothly  with  frequency  to 
500  cps,  where  a  resonance  peak  occurs.  'J?he  peak  may  be 
disregarded  as  being  caused  by  the  discontinuities  in  the 
correlation  functions. 

The  solid  line  in  Figures  22  and  23  is  the  theore- 
tical transfer  function  between  two  detectors.  Equation 
(92).  The  values  of  the  nuclear  parameters  used  in  SPAT 
are:   L^  =  1.61  cm^,  T  =  40  cm^,  /^  =  6  x  10"^  sec, 
/^  =  1  X  10*^  sec,  and  ^^=    0.923.  Only  /g  was  adjusted 
to  improve  the  agreement  between  the  theoretical  and 
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experimental  data.  The  other  parameters  were  obtained 
from  Reference  (29 ) . 

The  curves  in  Figures  22  and  23  were  not  determined 
by  using  the  least  squares  program.  It  is  apparent  that 
the  point  source  transfer  function  is  not  adequate  for  a 
light  water  and  uranium  medium.  The  large  size  of  T 
relative  to  L  makes  the  first  fli^t  collision  density- 
very  important  in  determining  the  thermal  neutron  distri- 
bution. In  any  future  work  using  Equation  (32)  with  a 
light  water  system,  the  expansion  coefficients  of  the 
source  (Equation  (15))  should  be  determined  for  a  dis- 
tributed source  (see  Reference  (40)  for  an  excellent 
example  of  a  source).  The  point  source  assumption  causes 
the  theoretical  transfer  function  to  predict  about  the 
right  attenuation  by  selecting  an  appropriate  Fermi  age, 
but  it  predicts  a  greater  phase  shift  than  actually 
occurs. 

Figures  24  through  30  present  the  results  from  the 
experiments  in  the  heavy  water  subcritical  assembly  with 
a  full  load  of  uranium. 

It  is  unfortunate  that  the  heavy  water  assembly  was 
constructed  next  to  a  4  foot  thick  concrete  wall.  The 
small  size  of  the  heavy  water  tank  (relative  to  L  and 
T)  combined  with  the  nearness  of  the  wall  and  floor  al- 
lows both  a  hi(^  neutron  leakage  rate  and  a  large  number 
of  reflected  neutrons  to  re-enter  the  system.  This  oc- 
curs even  though  the  tank  is  covered  with  cadmium.  The 
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reflection  peaks  in  Figure  24  show  the  strong  correlation 
between  the  original  source  neutrons  and  the  reflected 
neutrons.  The  fact  that  the  reflection  peaks  occur  at 
intervals  of  approximately  10  milliseconds  suggests  an 
explanation  for  the  heavy  water  data  presented  in  Figure 
16. 

The  Fermi  age  for  the  fast  neutrons  emitted  by  the 
source  is  large  enough  to  permit  almost  all  of  the  neu- 
trons to  leave  the  system.  Very  few  of  the  neutrons  ex- 
perience enough  collisions  to  become  thermalized  within 
the  systCTi.  The  escaped  neutrons  are  thermalized  and  re- 
flected back  to  the  tank  by  the  surroundings.  Those  neu- 
trons having  enough  energy  to  avoid  being  captured  by  the 
cadmium  re-enter  the  tank  and  complete  the  thermalizing 
process.  The  number  of  returning  neutrons  that  are  de- 
tected is  much  larger  than  the  number  of  original  neutrons 
detected.  Thus  the  strong  correlation  in  Figure  15  at  10 
milliseconds  in  a  reflection  peak.  Run  2  exhibits  only  a 
small  positive  correlation  at  T  =  0,  caused  by  the  detec- 
tors being  close  together.  Run  6,  where  the  detectors  are 
farther  apart,  exhibits  a  large  negative  correlation  at 
T  =  0,  indicating  that  the  variation  of  the  thermal  neu- 
tron density  is  quite  different  at  the  two  detector  loca- 
tions. 

As  uranium  is  added  to  the  system,  the  artificial 
source  produces  a  smaller  percentage  of  the  fast  neutrons, 
since  the  fast  neutron  source  is  spread  throughout  the 
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fuel.  This  spreading  of  the  source  reduces  the  leakage 
and  reflection  effect.  For  Run  40  (Figure  24),  the  cor- 
relation peak  at  T  =  0  is  stronger  than  the  reflection 
peak  at  10  milliseconds.  To  obtain  good  data  in  heavy- 
water  only  will  require  at  least  one  of  the  following 
items:  a.  a  larger  assembly,  b.  a  location  far  from  other 
structures,  and  c.  a  moderating  and  absorbing  structure 
around  the  tank  so  that  neutrons  are  not  returned  to  the 
tank. 

Figures  25  through  30  contain  the  data  used  to  deter- 
mine the  spatially  dependent  transfer  function  in  a  heavy 
water  subcritical  assembly.  The  portion  of  the  correlation 
curves  between  -3.5  and  +3.5  milliseconds  was  transformed. 
The  slight  difference  in  the  source  behavior  for  the  two 
runs  can  be  seen  in  the  auto  correlation  of  Detector  1, 
Figure  26.  The  power  spectra  data  appears  to  be  correct 
to  700  cps  (Figure  27),  with  the  discontinuity  effect  first 
appearing  at  250  cps.  The  positive  phase  shifts  shown  in 
Figure  28  are  caused  by  the  peaks  of  the  correlation  curves 
occurring  at  slightly  negative  delay  times  (Figure  25).  It 
is  not  coitpletely  clear  why  the  maximum  correlation  occurs 
at  this  time,  but  it  is  expected  that  the  pattern  and  timing 
of  the  neutron  reflections  and  thermalization  processes 
contribute  to  the  effect.  Taking  the  ratio  of  the  power 
spectra  (Equation  (91a))  cancels  the  reflection  effects  to 
some  extent,  as  seen  in  Figures  29  and  30. 

The  point  source  model  is  able  to  predict  reasonably 
well  the  attenuation  and  phase  shift  of  the  experimental 
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transfer  function  for  the  heavy  water  subcritical  assembly. 
The  theoretical  and  experimental  data  agree  fairly  well  in 
Figures  29  and  SO  for  frequencies  between  0  and  250  cps. 
The  resonance  peak  in  the  attenuation  curve  and  the  sudden 
increase  in  the  phase  shift  for  frequencies  greater  than 
250  cps  indicate  that  effects  other  than  simple  slowing 
down  and  diffusion  have  affected  the  correlations.  One 
of  these  effects  is  the  discontinuities  in  the  correlation 
curves.  Figures  25  and  26.  Another  possible  effect  is  the 
inhomogeneity  of  the  nuclear  system. 

Investigation  of  the  roots  of  the  characteristic  equa- 
tion for  the  theoretical  heavy  water  subcritical  assembly 
(Figure  13  and  Table  2)  shows  that  the  wave  length  of  the 
fundamental  mode  approaches  the  size  of  the  lattice  spacing 
in  the  actual  assembly  (12  cm)  as  the  frequency  increases 
(see  Figure  31  for  the  roots  at  0  and  2000  cps).  The 
higher  order  modes  have  wave  lengths  that  are  shorter  than 
the  lattice  spacing.  For  example,  the  wave  length  of  the 
fundamental  mode  at  50,  200,  and  400  cps  is  177,  50,  and 
30  cm.  The  wave  length  of  the  next  higher  mode  is  approxi- 
mately 8  cm  for  the  above  frequencies.  As  the  frequency 
is  increased,  the  att^iuation  of  the  fundamental  mode  in- 
creases faster  than  for  the  other  modes  (see  Figure  31), 
making  the  higher  modes  more  important  relative  to  the 
fundamental  mode.  Thus  the  nuclear  system  is  being  inves- 
tigated with  wave  lengths  of  the  same  order  of  magnitude  as 
the  lattice  spacing. 
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The  theory  contains  the  assumption  that  the  system 
is  homogeneous.  This  assumption  is  not  satisfied  in  the 
actual  experiment  at  high  frequencies.  The  sudden  change 
in  the  phase  shift  and  in  the  attenuation  (Figures  29  and 
30)  could  be  caused  by  the  lattice  acting  as  a  three  dimen- 
sional difraction  grating  (38.).  This  effect  is  not  seen 
in  the  light  water  subcritical  assembly  for  two  reasons. 
The  lattice  spacing  is  smaller  (4  cm)  and  the  higher  fre- 
quencies necessary  for  the  shorter  wave  lengths  are  not 
present  in  the  source  spectrum. 

The  values  of  the  nuclear  parameters  used  in  SPAT 
to  obtain  the  theoretical  curves  in  Figures  29  and  30  are: 
L^  =  92.9  cm^,  T  =  130  cm^,  ^g  =  0.95  milliseconds, 
f^-  0.1  milliseconds,  and  kgo  =  1.15.  L^,  T ,  and  k^ 
were  obtained  from  Reference  (39.).  /g  and^_,  were  ad- 
justed to  obtain  a  good  agreement  with  the  experimental 
data. 


CHAPTEl?  VI 
CONCUJSION 

The  experimental  results  clearly  show  a  spatial  effect 
that  must  be  included  in  any  theoretical  model  to  correctly 
describe  the  transfer  function  of  a  nuclear  system.  Vftiether 
the  model  used  to  develop  the  equations  in  Chapter  II  is 
adequate  is  not  definitely  known  at  this  time.  It  does 
appear  that  the  use  of  a  distributed  source,  based  on  the 
first  flight  collision  density  (see,  for  example,  T^eference 
(40) )  would  give  a  better  representation  of  the  light  water 
system,  'Aether  the  additional  labor  would  be  worthwhile 
for  other  systems  is  questionable. 

The  theoretical  results  show  that  the  lumped  para- 
meter model  is,  in  general,  not  adequate  for  describing 
the  transfer  function  of  a  nuclear  system.  However,  at 
a  specific  distance  from  the  source,  the  lumped  parameter 
model  does  predict  nearly  the  same  results  as  the  spatially 
dependent  transfer  function,  at  least  over  part  of  the 
frequency  range.  This  fact  may  be  used  in  selecting  the 
location  of  neutron  detectors. 

For  example,  the  theoretical  spatially  dependent  trans- 
fer function  equations  (3DTF)  developed  in  Chapter  II  can 
be  used  to  determine  the  locations  for  which  the  lumped 
parameter  transfer  function  is  reasonably  accurate  over  a 
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wide  range  of  k^f^  values.  If  the  detector  for  a  reacti- 
vity instrument  designed  to  utilize  Equation  (98)  is 
placed  at  one  of  the  locations  determined  by  SDTF,  the 
measured  impulse  response  function  can  be  easily  fitted 
using  the  simple  lumped  parameter  impulse  response  func- 
tion. Determining  the  reactivity  would  then  be  relatively 
easy.  If  the  detector  is  placed  elsewhere,  a  more  conqpli- 
cated  model  which  includes  the  spatial  effect  must  be  used. 

The  use  of  Equation  (91a),  which  comes  from  Equation 
(90),  to  determine  the  experimental  transfer  function  of 
a  system  has  several  advantages  over  other  statistical 
methods  (9,  10,  IS)   The  major  one  is  the  ability  to  cancel 
the  effects  of  the  source  and  the  effects  of  the  measuring 
system  from  the  final  result  by  taking  the  ratio  of  the 
cross  power  spectra.  Another  advantage  is  the  ability  to 
eliminate  the  effects  of  uncorrelated  noise.  The  main 
problem  in  using  the  cross  correlation  of  two  outputs  to 
measure  the  characteristics  of  a  system  comes  from  its 
accuracy.  Every  unwanted  effect,  such  as  the  reflection 
of  the  source  from  the  floor,  is  detected  and  added  to 
the  correlation  function.  The  correlation  function  cannot 
be  adequately  corrected  by  using  just  part  of  it  or  by 
shifting  it  vertically  to  make  the  function  equal  to  zero 
at  tN.  The  unwanted  effects  should  be  subtracted  from  the 
recorded  signals  before  correlating  them,  which  requires 
knowledge  of  the  undesired  signal.   It  would  be  better  to 
conduct  the  experiment  on  a  facility  specially  designed 
to  minimize  the  undesired  effects- 
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The  FM  tape  recorder,  which  was  acquired  after  the  ex- 
periments on  the  light  water  system  were  finished,  could 
be  used  in  future  work  with  that  system.  By  recording 
data  at  30  inches  per  second,  and  then  reproducing  the 
signal  at  1  7/8  inches  per  second  for  the  data  aquisition 
system,  a  time  expansion  of  sixteen  would  occur,  permitting 
the  fast  response  of  the  light  water  system  to  be  seen. 

The  use  of  a  pulse  detection  system  for  neutrons 
seems  to  be  much  superior  to  an  ion  chamber  and  current 
amplifier,  because  the  problems  of  60  cycle  per  second 
interference  are  greatly  reduced.   Other  noise  may  be  re- 
moved from  the  signal  by  use  of  a  discriminator  in  the 
linear  amplifier.  The  signal  supplied  to  the  count  rate 
circuit  is  essentially  noise  free.  Unfortunately,  further 
processing  of  the  signal  by  other  electronic  equipment 
soon  increases  the  noise  content  of  the  signal. 

Some  suggestions  for  further  work,  such  as  using  a 
distributed  source  to  find  the  source  expansion  coeffi- 
cients (Equation  (13)),  have  been  discussed  elsewhere 
when  applicable.  Another  possibility  for  further  work 
is  to  expand  the  theory  to  include  two  or  more  region 
systems,  using  the  approach  of  Foderaro  and  Garabedian 
(41).  Also,  the  transfer  function  (Equation  (46))  could 
be  transformed  on  the  U  variable,  rather  than  the  axial 
buckling  variable,  to  obtain  the  Fermi  age,  diffusion 
theory  impulse  response  function.  This  function  could 
be  used  directly  with  Equation  (98).  Further  study  of 
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the  neutron  reflection  phenomenon  and  its  effect  on  the 
correlation  function  could  prove  to  be  beneficial  and 
should  be  investigated  more  conpletely. 


APPENDICES 


APPENDIX  A 

FOTOIER  T!?ANgFORM  OF  A  PERIODIC 

DELTA  FaNCTION 

The  Fourier  transform  of  a  periodic  delta  function 
may  be  obtained  in  a  limiting  process,    starting  with 
periodic  pulses  of  width  b  and  height  a,    and  area  ab 
equal  to  one.     The  first  pulse  is  symetrical  about  t  =  0 
\^ile  the  rest  are  repeated  every     Tj  seconds.     In  in- 
tegral form,   the  Fourier  series  expemsion  is    (10.    12) 

00 

f(t)    =       ^      F(n)  exp(j    2TTnt/T-^)  (Al) 

n=-oo 

vAere  T]^  is  the  period  of  the  function  f  (t),  and  F(n), 
the  expansion  coefficients,  are 


^1/2 


(A2) 


F(n)  =  ^J     f(t)  exp(-j2  7rnt/T^)  dt 
-^1/2 


For  the  periodic  pulses,  F(n)  is  found  by 


j3/2 
F(n)  =  ^J   a  exp(-j2Trnt/T^)  dt         (A3) 
^-b/2 


X19 
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Since  f(t)  is  even.  Equation  (A3)  may  be  reduced  to 
the  cosine  transform  of  a,  which  is  (37) 


nTTb 

sm  -rrr 


^'"'  =   ^-^ 


(A4) 


In  the  limit,  as  a  goes  to  infinity,  b  goes  to  zero,  and 
ab  goes  to  1,  Equation  (A4)  becomes 


F(n)  =  ^ 
^1 


(A5) 


Using  this  result  in  Equation  (Al )  gives  for  the 
periodic  delta  function 


f(t)     =      ^    6(t.mT^) 


ra=-oo 
oo 


^^    exp(j27rnt/Ti)  (A5) 

n=-cD 


(31) 


APPENDIX  B 

COMBINATION  OF  TRIGONOTCTRY  FUNCTIONS, 

EQUATION  (50) 

The  function  in  Ecjuation    (50),    repeated  below,   may 
be  simplified  by  using  standard  trigonometry  functions 
as  given  in  Reference    (37,) . 

f(^)     =      y      cos^z  sin^  ^c  cos4mc^ 
m^-oD 

Using 

sin^^c     =     i-ico32^c  (B2) 

and 

cos  J  a  cos4cm^      =     ^  cos  ^(z  +  4cm)  +  ^  cos£  (z  -  4cm) 

in  Eguation    (Bl)  and  then  using  Equation    (B3)  again  in  the 
result,    f(^)  becomes 

CD 

f  ( J"  )     =      y  ^[2  cos  ^(z  +  4cm)   -I-  2  cos£  (z   -  4cm) 

m=-aD 

-  cos  ^  ( z  -f  4cm  +  2c )  -  cos  ^  ( z  +  4cm  -  2c ) 

-  cos^  (z  -  4cm  +  2c)  -  cos£  (z  -  4cm  -  2c )]     (B4) 

121 


122 

Consolidating  the  terms.  Equation  (B4)  is 


f(|)  =     J    Y.    f-l^^'cosICa  +  2cm) 

m=-oo 

00 

=  i  E  ^-l^"'exp[j£  (z  +  2cm)]      (B5) 


Either  form  may  be  used  in  Equation  (50) 


APPENDIX  C 
ROOTS  OF  THE  CHARACTERISTIC  EQUATION 

The  characteristic  equation  for  the  Fermi  age  theory, 
diffusion  theory  model  used  in  this  thesis  is  a  complex, 
transcendental  equation.  It  is  (see  Equation  (32)) 

DF  =  (H-  A^  ^.  l2^2^  jU^)  exp(j(Ji'^ 

with  the  eigenvalue,  or  buckling,  separated  into  the 
radial  buckling  and  the  axial  buckling.  Solving  Equation 
(CI)  for  the  axial  buckling  in  terms  of  the  frequency  and 
system  characteristics  allows  the  numerical  calculation  of 
Equation  (57),  the  inverted  transfer  function. 

Investigation  of  this  equation  for  the  steady  state 
case  has  been  done  by  others  {^).  Taking  the  roots  of 
Equation  (CI)  as 

^i  =   ±  jR^  =   ±j(A^  +  jB.)  (C2) 

with  1=0  corresponding  to  the  smallest  value  of  R. ,  it  is 
found,  for(J=0,  that  (JO.) 

^-i  =  ^ 


B..  =   -B, 
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0 

(C3) 
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or  the  roots  are  symetrical  about  the  imaginary  axis,  with 
the  locals  of  the  roots  crossing  the  imaginary  axis.  For 
CJfO,  the  locus  of  the  roots  is  not  symetrical.  As  an 
exaB5>le,  Figure  31  shows  a  plot  of  the  roots  for  a  fre- 
quency of  0  and  2000  cps.  The  parameters  used  in  Equation 
(CI)  were  the  same  as  those  used  in  calculating  the  trans- 
fer function  presented  in  Figure  1-^. 
Defining  ^f  as 

^\    =    ■   («i-^Jb,) 

=-[(A^-  b|)  -"-JZAjB.]  (C4) 

>2 

Equation  (C2)  may  now  be  written  in  terms  of  ^^. 


^i  =   -J 


(a^f  b^) 


J' 


[(af  -hbf)  - 


i 


(C5) 


>2 
Equation  (CI)  is  used  to  solve  for  Q^,   and  Equation  (C5) 

is  then  used  to  find  ^^. 

The  real  and  imaginary  parts  of  Equation  (CI)  may  be 

separated.  Substituting  Equation  (C4)  in  Equation  (CI), 

the  real  part  is 


(1+  L^B^  -  L^a^)  cos(Uje^  -  Tb.  )  -  (CJi 


-  lV)  siMui^  '  r 


v) 


exp(rB;: 


Ta.) 


k^[i^L;^rf(i->0)]    _ 


It  u^rf 


=  0 


(C5) 
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uio/x      ""v 
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and  the  imaginary  part  is 

2o2   t2 


1  +  u^r^ 


(C7) 


^ere  Q^    (Equation  (28))  has  been  simplified  by  consider- 
ing only  one  group  of  delayed  neutrons^  represented  by  '  j^. 

Equations  (C6  and  C7)  may  be  solved  numerically  for 
a^^  and  b , ,  using  standard  iterative  techniques. 

Nevton's  method,  as  applied  to  complex  equations,  is 
used  in  the  computer  code  described  in  Appendix  D,  Once 
an  estimate  of  the  value  of  a  root  is  obtained, 

t)PF 


A  ^2  _  .  __^i.^ (C8) 

^^         DF 


may  be  used  to  converge  on  the  correct  value  of  a  root. 
The  major  problems  are  insuring  that  R  is  found,  and  that 
no  roots  are  skipped.   For  most  systems,  only  the  first 
few  roots  will  be  needed  in  the  calculation  of  the  trans- 
fer function.  Equation  (57). 

A  good  first  guess  for  ^  may  be  obtained  one  of  two 
ways,  depending  on  the  values  of  the  parameters.  If  UA- 
and  TB^  are  greater  than  U(^  and  L  B^,  it  may  be  assumed 
that 


1+  l2b2+  l2^2  _^  JU/g     =     exp(L2B2+   l2^2^  jU/g 


)    (C9) 
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Using  the  definitions  (1,8) 


/t  =  fs^S. 


r 

Equation  (CI)  becomes 


and  the  solution  for  ^  is 


C 


m'L  KJ3 


^Mr-TT^^^n^^i^A 


(CIO) 


(Cll) 


(C12) 


Substituting  Equation  (CI 2)  in  Equation  (CI),  it  is 
seen  that  Equation  (C12)  is  an  exact  solution  if 

A     t 


and 


,/Jc  = 


(C13) 


mienUi       and  TB^  are  smaller  than  CJ  i     and  ITB^,    it 
/  "  s  n' 

may  be  assumed  that 


exp(r  ^^)  =  1  +  r^ 

Equation    (CI)  becomes 


(C14) 


f    [l+    l2b2+   jU4   .  Y^J3,  exp(.rB2   .  jUi^)]      =     0 


(CIS) 
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Equation  (C15)  is  a  second  order  equation  in  ^   and  may- 
be solved  using  the  quadratic  fornula. 

Both  Equations  (C12  and  C15)  are  available  in  the  com- 

/  2 
puter  code  for  calculating  the  first  guess  for  ^  .  As 

Ic   approaches  zero,  the  roots  of  the  non -multi plying 

medium  case  (Equation  (59))  may  be  used  for  the  first 

>  2  /  2 

guess  of  the  value  of  C  q»     Since  the  other  Q^   are  nearly 

periodic  in  the  b.  direction,  the  first  guess  for  them  may 


>  2 

be  based  on  Q^. 


APPENDIX  D 
COMPUTET^  CODES 

The  programs  used  to  calculate  the  spatially  dependent 
transfer  function  (Ecfuations  (57  and  63))  are  presented  in 
the  first  section  of  this  appendix.  The  second  section  con- 
tains two  programs  used  in  the  general  data  processing  code 
to  calculate  the  Fourier  transform  of  the  correlation  func- 
tion (Equation  (112)),  and  to  obtain  the  transfer  function 
(Equation  (91)). 

Transfer  Function 

Equations  (57  and  53)  are  evaluated  by  a  Fortran  II 
code  named  SPAT.  The  code  may  be  used  for  a  straight-for- 
ward calculation  of  the  transfer  function  as  a  function  of 
frequency  or  with  the  "General  Least  Squares  Program"  (24). 
If  two  distances  are  specified  in  the  input,  the  ratio  of 
the  transfer  functions  (Equation  (92))  is  calculated. 

The  main  program,  MAIN,  is  used  as  an  input  and  output 
program  when  the  least  squares  program  is  not  being  used  in 
that  capacity.  Subroutine  CALC;  has  an  argument  list  that  is 
compatible  with  the  least  squares  program.  All  input  is 
supplied  to  CALC  by  the  P  array,  which  is  defined  by  the' 
first  nineteen  instructions  in  CALC.  The  output  from  CAU: 
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is  controlled  by  P(4),  T^ich  is  defined  in  the  FINAL  sub- 
routine. Either  the  real  and  imaginary  part  of  the  trans- 
fer function,  or  the  phase  shift  and  the  attenuation  in  db, 
is  available.   Subroutine  GUE3T  calculates  a  normalizing 
factor  used  in  each  term  of  the  series  to  avoid  exceeding 
the  numerical  capacity  of  the  computer.  Subroutine  CRT 
finds  the  roots  of  DF  (Equation  (54)),  which  are  needed  for 
the  sum  over  i  in  Equation  (57).  The  guess  for  the  first 
root  is  made  by  GUESS,  using  the  approximate  equations  de- 
veloped in  Appendix  C.  Subroutine  TERM  evaluates  the  real 
and  imaginary  parts  of  each  term.  Subroutine  ERROR  calcu- 
lates the  magnitude  of  the  last  term  confuted,  relative  to 
the  sum  of  the  absolute  values  of  all  the  previous  terms. 
Both  the  real  and  imaginary  parts  must  be  less  than  a  speci- 
fied value  before  the  series  is  terminated.  The  first  fre- 
quency should  alvays  be  zero,  as  the  results  for  the  other 
frequencies  are  normalized  to  the  results  for  f  =  0. 
The  programs  are  listed  on  the  following  pages. 
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C      MAIN  PROGRAM  FOR  CALLING  CALC,  KYLSTRA 

DIMENSION   F(IOO),  P(40),  DO(AO) 
C      HAVE   P(4)  =   1.  FOR  THIS  PROGRAM 

CALL  BCDCON  (     50H   NF,   F(IOO),  P(40).  NC 
I  ,NF,  F,  P,  NO  ) 

C 

CL  =  8.6858896^ 
CT  =57.29578 
P(20)  =  0. 
P121)  =  0. 
1  CALL  SYMBOL   (  IN,  10) 
C      SYMBOL  IS  A  NON-FIXED  FORMAT  INPUT  SUBROUTINE 
C      SYMBOL  AND  BCDCON  MAY  BE  REPLACED  BY  REGULAR  INPUT 
E  =  P(3) 
CK  =  P(A) 
NP  =  P(8) 
FN  =  P(  16) 
CFAC  =  P(17) 
EF  =  P( 18) 
DFAC  =  P(19) 
5  WRITEOUTPUTTAPE  6,  1001 
1001  FORMAT  (IHl/  37X  28HRESULTS  FOR  MULTIPING  MEDIUM 
1   18H  TRANSFER  FUNCTION     ) 
WRITEOUTPUTTAPE  6.  1003, P(l),  P(2),  P( 7  )  , P (9 ) , P ( 10 ) 
1003  FORMAT  (  IHO   14X      26HVALUES  OF  PARAMETERS  ARE. 
1   9HRADIUS  = 

1  F8.1,   18H  CM.  AND  LENGTH  =   F6.1, 

2  15H  CM.  FOR  TANK,    /   15X  9HRADIUS  =   F6.1, 

3  19H  CM.  AND  LENGTHS  =   2F6.1, 

4  28H  CM.  FOR  DETECTORS  2  AND  3.   /   ) 
WRITEOUTPUTTAPE  6,  1005,  P(15).  Pill),  P112),  P(14) 

1,  P(13),  P15),  P(6) 
1005  FORMAT  {  IH   14X  20HMULT I  PL ING  FACTOR  =   F7.4, 

1  21H,  DIFFUSION  LENGTH  =   F6.1,  18H  CM.,  FERMI  AGE  = 

2  F6.1,  6H  CM2.   /  15X  29HTHERMAL  AND  FAST  LIFETIMES  = 
32F9.6,  27H  SECS.,  AND  BETA  AND  TAU  =   F9.6,  F8.2, 
4//7X  AHFREQ  12X  AHREAL  12X  4HIMAG  13X  3HAMP  14X  2HD8 

5  12X  5HTHETA  IIX  5HTHETA  6X  5HTERMS    ) 
CO   AO  I  =  I,  NF 

W  =  F(I) 

CALL  CALC  I  W,  XI,  X2,  A,  P.  DO,  XK  ) 
Al  =   SQRTF  (  A»A  +  XK»XK  ) 
A2  =   CL  ♦  LOGF  (Al) 
THE  =  CT  •  ATANF  (  XK  /  A  ) 
TH2  =  180.  -  THE 

WRITEOUTPUTTAPE  6,  1004, W, A , XK, Al , A2 ,    THE, TH2, DDl 1 ) 
1004  FORMAT  {  IH  7E16.7  ,  F7.0  ) 
40  CONTINUE 
GOTO  1 
END 
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C      SPATIALLY  DEPENDENT  TRANSFER  FUNCTION,  KYLSTRA 

C      FOR  A  MULTIPLING  MEDIUM,  POINT  SOURCE 

C      FOR   A   NON  -  MULTIPLING  MEDIUM,   POINT  SOURCE 

C      P  VECTOR  DEFINED  BY  FIRST  19  INSTRUCTIONS 

C      XK  CONTAINS  UNWANT  PART  OF  H,   REAL  OR  IMAG 

C      OR  DB  AND  THETA 

C 

SUBROUTINE  CALC  (  F,  XI,  X2 ,  H,  P,  DO,  XK  ) 

DIMENSION   P(40),  BE(20),  BZ(20)  ,  00(40) 

R    =   P(I) 

C    =   P(2) 

E    =   P(3) 

NCK  =   P(4) 

B    =   P(5) 

Tl   =   P(6) 

RO   =   P(7) 

NP   =   P(8) 

12      =   P19) 

Z3   =   P(10) 

DL   =   Pill) 

T    =   P(12) 

SL   =   P(I3) 

TL   =   P(U) 

RKF  ^       PI  15) 

FNN  =   P( 16) 

CFAC  =  P117) 

EE   =   PdS) 

DFAC=   P{19) 
C      CFAC  FOR  DELTA  B,  OFAC  FOR  DELTA  R  AND  S 

IF(P{20))  550,  550,  555 
550  CONTINUE 
C      20  BESSEL  FUNCTION  VALUES 

BEd)   =   .269465 

BE(2)   «   .115804 

BE13)   =   .073712 

BE(4)   =   .054056 

BE(5)   =   .042642 

BE(6)   =   .035231 

BE17)   =   .030033 

BE(8)   =   .026169 

BE{9)   =   .023165 

BEtlO)  =   .020794 

BE( 11)  =    .018851290 

BEI12)  =    .017241414 

BE{13)  =    .015901210 

BEI14)  =    .014737960 

BE{  15)  =    .013735840 

BE(16)  =    .012859560 

BE(17)  =    .012100000 

BEI18)  =    .011406240 

BE( 19)  =    .010816000 
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BE[20) 

.010261690 

BZ(1) 

=   2.4048 

BZ(2) 

=   5.5201 

BZ13) 

=   8.6537 

BZ(4) 

=  11.7915 

BZ(5) 

=  14.9309 

BZ{6) 

=  18.0711 

BZ(7) 

=  21.2116 

BZ(8) 

=  24.3525 

BZ(9) 

=  27.4935 

BZI 10) 

=  30.6346 

BZl 11) 

=  33.7758 

BZ(12) 

=  36.9171 

BZl 13) 

=  40.0584 

BZl  K) 

=  43.1998 

BZl  15) 

=  46.3412 

BZl 16) 

=  49.4826 

BZl  17) 

=  52.6241 

BZl  18) 

=  55.7655 

BZl  19) 

=  58.9070 

BZ120) 

=  62.0485 

AN  =  BN 

RC2  =  ! 

5C2 

RC4  =  SC4 

PI2  =  6.28318531 

PJ20)  = 

-    1. 

NX  =  P121) 

NY  =  P121) 

NPl  =  0 

555  CONTINUE 

DL2   = 

DL     •»  2 

W 

PI2  *  F 

ws  = 

W  •  SL 

WT   = 

W  »  TL 

S3  =  0. 

R3  =  0. 

S2  =  0. 

R2  =  0. 

56  =  0. 

R6  =  0. 

AR3  =  0. 

AS3  =  0. 

AR2  =  0, 

AS2  =  0. 

CL  =  2.»C 

DB 

1.  -   B 

TIW  = 

W  •  Tl 

T1W2  = 

TIW  •♦  2 

RKD   = 

RKF    /  (  1.  +  TIW2) 

RK   = 

RKO  »  I  1.  +  T1W2  •  DB) 

SK   = 

RKD  ♦  TIW  •  B 
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COU  =  -1. 

.IF  (  F  )  600,  600,  625 
600  IF(  NY  )  605,  605,  300 
605  NY  =  1 

IFIRKF-  .05)  510,  510,  500 
500  CONTINUE 

B2S  =  (  2.4  /  R)  »»2 

CALL  GUEST  {  B2S,  DL2,  T,RKF,  RC2,  SC2  ) 

CRR  =  Z3  •  SQRTFC  RC2) 
;      CONSTANT  FACTOR  IN  EACH  TERM 

GO  TO  300 
510  CRR  =  Z3  /  OL 

GO  TO  300 
625  NY  =  0 

300  CONTINUE 

DO  70  N  =  1,  20 

B2N  =   (  BZ(N)/  R  )  ♦•  2 

DLBN  =   DL2  ♦  B2N 

TBI  =  T»B2N 

AA  =  1.    . 

301  CONTINUE 

IF(  RO  )   3,  3,  2 

2  AA  =       BJOC  BZ(N)»  RO  ) 

:      BJO  CALCULATES  THE  JO  BESSEL  FUNCTION 

3  CONTINUE 

COl  =  1.  +  DLBN 
C02  =  AA  /  BE{N) 
NR  =  0 

4  CONTINUE 
EF  =  EE 

CALL  CRT  (C01,DL2,TB1,T,WS,WT,RK,SK,NP1,NR,AD2,B02, 
IAN,BN,EF,CFAC,DFAC) 
IFIEF  -  9.5)  302,  302,  2001 
2001  H  =  0. 
XK  =  0. 
GO  TO  101 

302  CONTINUE 

RCl  =   T  ♦  A02  -  TBI 
SCI  =   T  «  BD2 

CALL  CEX  {RCl,  SCI,  RC2,  SC2) 
IFIRK)  303,413,303 
A13  RC3  =  1. 
SC3  =  0. 
GO  TO  415 

303  CONTINUE 

RC3  =  DL2  ♦  T  •  {  COl  -  DL2»  AD2 ) 
SC3  =   T  *  (  WT  -  DL2  ♦  BD2  ) 
415  CONTINUE 

CALL  CMUL  (AN,  BN,  RC3,  SC3,  RC4,  SC4  ) 
CALL  COIV  (  RC2,  SC2,  RC4,  SC4,  RC5,  SC5) 
RN  =  RC5  *  C02 
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SN  =  SC5  •  C02 

THESE  ARE  CONSTANTS  FOR  M  INDEX 
M  =  0 
FM  =  0. 
NCT  =  0 
NCTl  =  0 
GT  =  0. 
505  CONTINUE 

CALL  TERM  I Z 3, GT, M, AN, BN.RN, SN, DR , DI ,CRR) 

S3  =  S3  ♦  DI 

R3  =  R3  +  DR 

CALL  ERROR!  DR  ,  DI  ,  AR3,  AS3,  RAR3,  RAI3  ,  NR) 

IF(Z2)  405,A05,5 

5  CONTINUE 

CALL  TERM  ( Z2,GT,M, AN, BN,RN, SN,DR,DI ,CRR) 
R2  =  R2  +  DR 
S2  =  32  ♦  DI 

CALL  ERRORl  DR  ,  DI  ,  AR2,  AS2,  RAR2,  RAI2  ,  NR) 
405  CONTINUE 

COU  =  COU  +  I. 

6  IF{  COU)  30,  30,  9 
9  CONTINUE 

IF(  NCT  -  100)  55,  55,  50 
50  NCT  =  2 

NCTl  =  NCTl  +  1 

WRITEOUTPUTTAPE  6,  1001 , NCT 1 , NR, N, R3 , S3 ,R2, S2 , F , 
1  RAR3,RAI3,RAR2,RAI2 
1001  FORMAT!  41H0  VALUES  OF  NCTl,  NR,  N,  R3,  S3,  R2 ,  S2, 

1  31HF,  RAR3,  RAI3,  RAR2,  RAI2.  ARE    /3I5,   4620. 8, 

2  /  5E20.8   ) 

IF  INCTl  -  5)   55,  55,   75 

55  IF  I  RAR3  -  E)  56,56,  30 

56  IF  (  RAI3  -  E)  57,57,  30 

57  IFl  Z2)   65,  65,  60 

60  IF  (  RAR2  -  E)  61,  61,  30 

61  IF  {  RAI2  -  E)  65,  65,  30 
30  FM  =  FM  +  1. 

M   =   M  +  1 
NCT  =  NCT  +  1 
GT  =  CL  •  FM 
:      SKIP  NEXT  M  IF 

IF  l(  GT-  2.»Z3)»AN  -  10.  )   505,  65,  65 
65  IF{RK)   201,  201,  250 
201  IF  (FM)   75,  75,  70 
250  IF{  FM    )  260,  260,  4 

260  IF  (  NP  )   262,  262,  261 

:      NP  =  1,  CRT  PRINTS  OUT  ROOTS 

261  CALL  CRT  { COl ,DL2, TBI , T, WS . WT ,RK,SK, NP, NR, AD2 , BD2 , 
1AN,BN,EF,CFAC,0FAC) 

262  CONTINUE 

IFtNR  -  1  )  75,  75,  70 
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70  CONTINUE 
C      END  OF  N  DO  LOOP 
75  CONTINUE 

IF(  Z2  )  80,  80,   85 
80  CONTINUE 

RR5  =   COSFI  WS) 
SR5  =  -SINF{  WS) 

CALL  CMUL  (  RR5,  SR5,  R3,  S3,  R6,  S6  ) 
GO  TO   89 
85  CONTINUE 

CALL  CDIV  (  R3,  S3,  R2,  S2,   R6,  S6  ) 
89  CONTINUE 
C      PICK  FIRST  F  =  0  FOR  NORMALIZING 
C      AND  COMPONENT  DESIRED  FOR  LEAST  SQUARES 
CALL  FINAL  (  F,  R6,  S6,  NCK,  H,  XK) 
101  CONTINUE 

DD{  1)  =  COU   +  I. 

RETURN 

END 


SETS  UP  CRR  FOR  NORMALIZING  EACH  TERM,  KYLSTRA 
SUBROUTINE  GUESTI  B2S,  DL2,  T,  RK,  RC2,  SC2  ) 
TBI  =  T»B2S 


DT  =  1.  +  DL2»B2S 

A  =  DL2  •  T 

B  =  I  T»DT  +  DL2  )  /  . 

(  A  »  2,) 

C  =  (  DT  -  RK*EXPF  {  - 

-TBI  ))/A 

Q  =  B»B  -  C 

IFl  Q  )  10,  20,  20 

10 

RC2  =   B 
GO  TO  50 

20 

RC2  =  B  -  SQRTF  I  Q  ) 

50 

SC2  =  0. 

RETURN 
END 

CALCULATES  EACH  TERM,   KYLSTRA 

SUBROUTINE  TERM  (Z ,GT, M, AN, BN,RN , SN, DR, DI , CRR) 

GP  =   Z  +  GT 

ZFP  =  GP»AN  -  CRR 

CALL  CEX  (-  ZFP    ,  -  GP»BN,  Rl,  SI) 

IF  {  M  )  20,  20,  5 

GP  =   GT  -  Z 

ZFP  =  GP»AN  -  CRR 

CALL  CEX  {-  ZFP    ,  -  GP*BN,  R2,  S2) 

SIG  =  (-l.)»»  M 

Rl  =  SIG  *  (  Rl  ♦  R2  ) 
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SI  =  SIG  »  (  SI  +  S2  ) 

20  CONTINUE 

CALL  CMUL  I  RN,  SN,  Rlt  SI,  OR,  01) 

RETURN 
END 


CALCULATES  ABS.  ERROR,   KYLSTRA 

SUBROUTINE  ERROR  (  OR,  DI,  R,  S,  RAR,  RAI,  NR) 

DAR  =  ABSF(DR) 

DAS  =  ABSFIDI) 

R  =  R  ♦  DAR 

S  =  S  +  DAS 

IF{  XMODF  {NR,2))  50,  20,  50 
20  RAR  =  1. 

FORCES  TERMINATION  ON  ODD  ROOT 

GO  TO  100 

50  RAR  =  DAR  /  R 

100  RAI  =  DAS  /  S 

RETURN 

END 


C      SUBROUTINE  CRT  FOR  COMPLEX  ROOTS  OF  D.F.,  KYLSTRA 

SUBROUTINE   CRT  (  DBl,  DL,  TBI,  T,  WS,  WT,  RK,  SK,  NP 
1,  NR,  A,  B,  AA,  BB,  E,  F,  F2) 
C      DBl  =  1+L2*32,DL  =  L2,  TBI  =  T»B2,  T  =  T,  WS  =  W*  SL, 
C      R;<  =  K»l  1  +  W2»T12»(1-B))/1  +  W2»T12,  SK=K»W»Tl»B/ l  +  W2»T12 
C      NP  =  1,  PRINT  ROOTS,  =0,  DONOT,  NR  =  ROOTS  FOUND 
C      A  +  J  B  =  ROOTS,  AA  +  J  BB  =  SQRT  (A  +  J  B  ) 
C      F  =   FACTOR  ON  TWO  PI,  F2  =  FACTOR  ON  DR  AND  DI 
C      WT  =  W»TL,E  =  ERROR  ALLOWED 

DIMENSION  AR{  IGO ) , BR ( 100 ) , AAR (100) ,BBR ( 100) ,COUT ( 100) 

KC   =  0 

F3  =  F 

Fl  =  F2 

IF(  NP)   20,  20,   2 
C      NP  =  1,  PRINT  OUT  ROOTS 
2  CONTINUE 

RC3  =  SC3 

RC4  =  SC4 

YP  =  ZP 

DR  =  DS 

WRITEOUTPUTTAPE  6,  1001 , DBl , DL, TBI , T ,WS ,WT ,RK, SK 
1001  FORMAT  (  IHO   20X34HVALUES  OF  ROOTS  FOR   THE  FOLLOWING 
1  lOH  CONSTANTS 

1  //   9X  3HLB1  i2X  2HDL  UX  3HTB1  13X  IHT  12X  2HWS  12X 

2  2HWT  12X  2HRK 

2  12X  2HSK  /    AX   8E14.6,  //   UX   AHROOT  13X  IHA 

3  15X  IHB  14X 
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32HAA   14X   2HBB     12X   5HC0UNT  /  ) 
DO  5  I  =  1,'  NR 

WRITEOUTPUTTAPE  6,  1002,  I,  AR{I),  BR(I),  AARd), 
I  BBR(I) ,  COUT(I) 
1002  FORMAT  (  IH   9X  15,  5X  5E16.7  ) 
5  CONTINUE 
GO   TO   600 
20  CONTINUE 

IF{NR)   25,  25,    100 
C      IF  1  ROOT  ALREADY  FOUND,  PROCEED  FROM  IT 
25  CONTINUE 

SKT  =  6.28  •  F3  /  T 
IF{RK  -  .05)  30,  30,  50 
C      IF  -,  USE  DIFFUSION  APPROX,  OTHERWISE,  GUESS  ROUTINE 
30  Rl  =      DBl  /  OL 

51  =      WT  /  DL 
C      EQUATION  58 

IF  (RK)   500,  500,  200 
C      RK  =  0,  GO  TO  SORT  PART,  OTHERWISE,  TO  CONVERGENT  PART 
50  CONTINUE 

D32  =  DBl 

DL2  =  DL 

TB2  =  TBI 

T2  =  T 

WS2  =  WS 

WT2  =  WT 

RK2  =  RK 

CALL  GUESSJ  DB2,  DL2,  TB2,  T2 ,  WS2,  WT2 ,  RK2,  Rl,  SI) 

GO  TO   200 
100  CONTINUE 
C      ESTIMATE  NEXT  ROOT  FROM   LAST 

52  =   BRINR) 

50  =  BR(1) 

51  =  -  S2  +  SO 
IFlNR-2)  90,  90,  95 

90  Rl  =  AR( 1) 
GO  TO  96 

95  Rl  =  AR{NR-1) 

96  CONTINUE 

IF{   XMODF  (  NR,  2  ))   105,  200,   105 
105  CONTINUE 

SI  =  SI  +  SIGNFISKT,  SO) 
C      END  OF  GUESSING  SECTION,  START  CONVERGENT  SECTION 

200  CONTINUE 
COU   =  0. 

201  IDY   =   1 
R  =   Rl 

S  =   SI 

IF(  Fl  )  210,  207,210 
207  Fl  =  1. 
210  CONTINUE 
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RCl  =  DBl  -  DL  *  R 
SCI  =  WT   -  DL  »  S 
RC2  =  TBI  -  T   »  R 
SC2  =  WS   -  T   •  S 
CALL  CEX  (  RC2,  SC2,  RC3,  SC3  ) 
CALL  CMUL(  RCl,  SCI,  RC3,  SC3,  RC4,  SC4  ) 
Yl=   RC4  -  RK 
Zl=   SCA  +  SK 
COU  =  COU  +1. 
NCT  =  NCT  +  1 
RC5  =   DL  +   T»  RCl 
SC5  =   T  »  SCI 

CALL  CMUL  1  RC3,  SC3,  RC5,  SC5,  YP,  ZP  ) 
YP  AND  ZP  ARE   DERIVATIVES   OF   Y  AND  Z 
CALL   CDIV  (Yl.Zl.  YP,  ZP,  DR ,  DS  ) 
Rl  =  R   +  DR  ♦Fl 
SI  =  S   +  DS  »F1 
250  CONTINUE 

RE  =   ABSFl  OR/  R  ) 
SE  =   ABSFt  DS/  S  ) 
IF{  RE  -  E  )  255,  255,  260 
255  IF(  SE  -  E  )  270,  270,  260 
260  R  =  Rl 
S  =  SI 

IFtCOU  -   99.)  210,  210,  268 
RETURN  FOR  NEXT  ITERATION 
268  KC  =   1 

GO  TO   2 
270  CONTINUE 
500  A  =  Rl 
B  =  SI 
NR  =  NR  +  1 
AR(NR)  =  Rl 
BR[NR)  =  SI 

CALL  CSRT  (  Rl,  SI,  RCl,  SCI) 
AA  =  RCl 
BB  =  SCI 
AAR(NR)  =  RCl 
BBRINR)  =  SCI 
COUT(NR)  =  COU 
IF(  NR  -100)   610,510,  510 
510  KC  =  1 

GO  TO  2 
600  IF  (  KC)    610,   610,  605 
605  CONTINUE 

WRITEOUTPUT  TAPE  6,  1005,  KC,  NR ,  COU 
1005  FORMAT  (  IHO  15X  6H  KC  =   15,  7H.  NR  =   15, 
1   8H.   COU  =    F7.1,  /  15X  IIHSKIP  ONE      ) 
E  =  10. 
610  RETURN 
END 
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C      SUBROUTINE  GUESS,  KYLSTRA 

C      USED  BY  CRT  FOR  FIRST  ROOT  OF  DL 

SUBROUTINE  GUESSl  DBl,  DL,  TBI,  T,  WS,  WT,  RK,  Rl,  SI  ) 

IF(  WS  -  WT)  1,  2,  2 
C      WS  GREATER  THAN  WT,  USE  EXP  SCHEME 

1  NW  =  1 
WFAC  =  WT 
GO  TO  5 

2  NW  =  -1 
WFAC  =  WS 

5  CONTINUE 

DT  =  DBl  -  1. 
IF{  TBI  -  DT)  7,  8,  8 
C      TBI  GREATER  THAN  DT,  USE  EXP  SCHEME 

7  NC  =  I 
TFAC  =  DT 
GO  TO  10 

8  NC  =  -1 
TFAC  =  TBI 

10  CONTINUE 

IF  (  WFAC  -  TFAC  )     12,  13,  13 
C      GO  TO  LARGEST 

12  NU  =  NC 
GO  TO  15 

13  NU  =  NW 
15  CONTINUE 

IF  (  NU)  20,  20,  50 
C      NU  +,  USE  POLY,  NU  -,  USE  EXP 
20  TIG  =  DL  +  T 

Rl  =  (  -LOGF{RK)  +  DT  +  TBI  )  /  TIG 
SI  =  (WS  +  WT)/TIG 
GO  TO  100 
50  CONTINUE 

TLF  =  2.  *T  *DL 
PI  =  T»DB1  +  DL 
Ql  =  WT»  T 
P2  =  Q2 


CALL  CEX  {  -TBI,  -WS,  P2,  Q2) 

P3  =  DBl  -  RK*P2 

03  =  WT  -  RK»Q2 

P4  =  Q4 

CALL  CMUL  {  PI,  Ql,  PI,  Ql,  P^," 

Q4) 

P5  =  P4  -  2.  »  TLF  ♦  P3 

Q5  =  Q4  -  2.  »  TLF  •  Q3 

CALL  CSRT  (  P5,  Q5,  P3,  Q3) 

Rl  =  (PI  -  P3)  /  TLF 

SI  =  (Ql  -  Q3)  /  TLF 

100 

CONTINUE 

RETURN 

END 
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C  FINAL    NORMALIZES    AND    SETS    UP    L.S.    COMPONENT.       KYLSTRA 

SUBROUTINE 'FINAL(F,    R5,     S5,    NC,    T,     X) 
DBLF(X,Y}  =  4.3429AA8  ♦  LOGF  (  X«X  +  Y»Y) 
R6  =  R5 
S6    =    S5 
NCK  =  NC 

IF     (    F    )    520,     520,     525 
520    IF     (    NX    )    528,    528,     530 
528    RNF    =    R6 

AMPN    =    DBLF    (    R6,    0.) 
Na  =  I 
GO    TO    530 
525    NX    =    0 
530    CONTINUE 
C  RETURN    COMPONENT    DESIRED    FOR    LEAST    SQUARES    FIT 

GO    TO    1     535,    535,    550,     550    ),    NCK 
C  CK    =    1,    REAL,    =    2,     IMAG,    =    3.    DB ,    =    A,     THETA 

535    R6    =    R6    /    RNF 
S6    =    S6    /    RNF 
GO    TO    I     540,    545    ),    NCK 
540    H    =    R6 
XK    =    S6 
GO    TC    101 
545    H    =    S6 
XK    =    R6 
GO    TO    101 
550    AMP    =    DBLF(R6,S6)       -    AMPN 

TH    =    57.29578    *    ATAiMF    (     S6/R6) 
GO    TO    (     540,    545, 555, 560), NCK     . 
555    H    =       AMP 
XK    =    TH 
GO    TO    101 
560    XK    =    AMP 

H    =    TH 
101    CONTINUE 
T    =    H 
X    =    XK 
RETURN 
END 
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C      SUBROUTINE  CMUL,  KYLSTRA 
C      SOLVES  R3  +  J  S3  =   (  Rl  + 

SUBROUTINE  CMUL  (  Rl,  SI, 

R3  =   Rl  •  R2  -  SI  »  S2 

S3  =   Rl  »  S2  +  R2  •  SI 
RETURN 
END 


J  SI  )  »  (  R2  +  J  S2  ) 
R2,  S2,  R3,  S3) 


c 

SUBROUTINE  CDIV,  KYLSTRA 

c 

SOLVES  R3  +  J  S3  =  (  Rl  +  . 

J  SI  )  /  (  R2  +  J 

SUBROUTINE   CDIV  I    Rl,  SI, 

R2,  S2,  R3,S3  ) 

DEM  =  R2»R2  +  S2  »  S2 

R3   =  {Rl  •  R2  +  SI  »  S2  ) 

/  DEM 

S3   =  (R2  «  SI  -  Rl  *  S2  ) 

/  DEM 

RETURN 

END 

c 

SUBROUTINE  CSRT  ,   KYLSTRA 

c 

SOLVES  R2  +  J  S2  =  SQRTF  I 

Rl  +  J  SI  i 

SUBROUTINE  CSRT(  Rl,  SI.  R2,  S2) 

R3  =  Rl 

S3  =  SI 

COl  =  SQRTF(  R3»R3  +  S3*S3 

) 

R2  =  SQRTF   ({  COl  +  R3  )  •.S  ) 

SA  =  SQRTFi (  COl  -  R3J  *,5 

) 

100  S2  =  SIGNF  {  S4,  S3) 

RETURN 

END 

c 

SUBROUTINE  CEX,  KYLSTRA 

c 

SOLVES  R2  +  J  S2  =  EXPF  {  Rl  ♦  J  SI  ) 

SUBROUTINE   CEX  (  Rl,  SI,  R2,  S2  ) 

SA   =  EXPF{  Rl) 

R2  =  SA  *  COSFl  SI) 

S2  =  SA  »  SINF(  SI) 

RETURN 

END 

S2  } 
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Fourier  Transform 

A  Fortran  II  subroutine  named  TRNS  evaluates  Equation 
(112).  It  is  called  by  the  general  data  processing  program 
(31 )  to  calculated  the  Fourier  transform  of  the  correlation 
functions.  The  majority  of  the  input  and  output  from  THNS 
is  through  the  common  storage;  however,  the  frequencies  at 
vhich  the  transform  is  calculated  may  be  read  in.  The  fol- 
lowing pages  contain  a  listing  of  TRNS  and  NXTPR,  an  output 
subroutine  used  to  combine  the  results  of  successive  runs 
for  H23(CJ)  (Equation  (91)). 
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:      TRNS,  LINEAR  INTEGRATION  OF  FOURIER  TRANSFORM,  KYLSTRA 

SUBROUTINE  TRNS 
:      COKMCN  FOR  VARIABLES  SAME  IN  DIFFERENT  MZ  CASES 

COMMON  X,H,G,NCT,ID,  I  D2 , FNT  ,NRD, NDS , M, MZ ,MM , MZ2 . MPL 

COMMON  NFILE,WMAX,NDEX,NFREQ,GT 

COMMON  VCTYI,VCTYK,VCTAC 

COMMON  JVAjJVB 

DIMENSION  XI 600 3 ,Hl60n ,G{600) 

DIMENSION  VCTVI(40),VCTYK{40) ,VCTAC(40) 

DIMENSION  GT(  18) 

DIMENSION  JVA{40), JVB{40) 
;      COMMON  FOR  VARIABLES  THAT  DIFFER  FROM  MZ  TO  MZ 

COMMON  Y,AC 

COMMON  TRR.TRI 

DIMENSION  YI750,A),  AC(300,4,4),  TRR1300,4,4) 

DIMENSION   WKW1300),  TRI{300,4,4) 

EQUIVALENCE  (  X,  WWW) 

TA  =  NRD  »  MZ 

TA2  =  TA  /  15000. 

FT  =  MM 

TAM  =  TA2  »  FT 

READINPUTTAPE  5,  5001,  NPRT,  NCAL,  NTRY,  NWIN 
5001  FORMAT  (  515  ) 

WWW( I)  =  0. 

NCTl  =  1 

IF  (  NPRT)   1980,1990,  1980 
1980  CONTINUE 

NCT1=  NPRT  +  1 

READ  INPUTTAPE  5, 1981 , ( WWW( K) ,  K  =  2,  NCTl) 


1981 

FORMAT  (   5E12.1  ) 

1990 

CONTINUE 

IF{  NCAL  )  1997,  1997, 

1991 

1997 

IF(  NPRT  )  1999,  1998, 

1999 

1998 

NCAL  =  M 

1991 

FI  =  I. 

NBB  =  NCTl  +  1 

NCTl  =  NCTl  +  NCAL 

DW  =  .5  /{  TA2  »(FLOATF(  NCAL  -  1  ))) 

DO  1995  KK  =  NBB,  NCTl 

WWWIKK)  =  DW  »  FI 
1995  FI  =  FI  +  1. 
1999  NMT  =  NCTl 

WRITEOJTPUTTAPE  6,  1940,     (WWWtll),  II  =  l.NCTi: 
1940  FORMAT  (  IHO  15X  34HTHE  SPECTRAL  DENSITY  COEFFICIENTS 
145HT0  FOLLOW  ARE  EVALUATED  AT  THESE  FREQUENCIES 
1  //  (  4X  5E20.8  )) 

DO  3000   L  =  1,  NDEX 

J  =   JVAIL) 

K  =   JVB(L) 

IT  =  K  -  J 
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FIT  =  IT 

T  =  FIT  /  15C00. 
C      FIRST  DELAY  POINT  IS  T  FOR  CROSSCORRELAT ION 
C      FIRST  DELAY  POINT  IS  ZERO  FOR  AUTOCORRELATION 

TLR  =   AC(M,J,K)  -  AC{MM,J,K) 

TLI  =   AC(M,J,K) 
C      LAST  POSITIVE  TAU  TERMS 

TCR  =  AC(I,J,K)  -   AC(2,J,KJ 

IFIIT)  2001,  2005,  200L 
2001  TCR  =  TCR  +  AC{1,K,J)  -  AC{2,K,J) 
C      CENTER  TERMS 
2003  TFI  =  AC(M,;<,J) 

TFR  =   AC;M,K,J)  -  AC(MM,K,J) 
C      LAST  NEGATIVE  TAU  TERMS 

2005  CONTINUE 
C      START  ZERO  FREQUENCY  INTEGRATION 

RHO  =   ACIM,J,K)  *  .5 

DO   2010  KK  =  1,  MM 

2010  RHO  =   RHO  +  AC[KK,J,K) 
IF{IT)  2011,  2014,  2011 

2014  RHO  =  RHO  -  .5  »  AC(1,J,K) 
GO  TO  2013 

2011  RHO  =  RHO  +  AC(M,K,J)  ♦  .5 
DO  2012   KK  =  2, MM 

2012  RHO  =  RHO  +  ACIKK,  K,J) 

2013  CONTINUE 

STRl  =  RHO  »  TA2 

H(L)  =  STRl 

TRR{1,J,K)  =  1. 

TRI  {1,J,K)  =  0. 
C      HAVE  FREQ  EQUAL  ZERO  RESULT 

C      CALCULATE  AND  STORE  CHANGE  OF  SLOPE  IN  AC  MATRIX 
C      ACID  =    AC(3)  -  2AC(2)  •»•  AC(1) 

PHI  =  AC(1, J,K) 

PH2  =  AC(2,J,K) 

DPI   =   PH2  -  PHI 

DO  2020   KK  =  3,  M 

KF  =  KK  -  2 

PH3  =  AC(KK,  J,  K) 

DP2  =  PH3  -  PH2 

ACIKF,  J,K)  =   DPI  -  DP2 

DPI  =  DP2 

PH2  =   PH3 
2020  CONTINUE 

IFI  IT)   2025,  2050,  2025 
2025  PHI  =   AC(  1,K,J) 

PH2  =   AC{2,K,J) 

DPI  =   PH2  -  PHI 

00   2040   KK  =  3,  M 

KF  =  KK  -  2 

PH3  =   ACIKK,  K,  J  J 
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DP2  =   PH3  -  PH2 

AC{KF,K,J)  =  DPI  - 

DP2 

DPI  =  DP2 

PH2  =  PH3 

2040 

CONTINUE 

2050 

CONTINUE 

THIS  LOOP  CONTROLS 

FREQ 

DO   2500  KK  =  2,  NMT 

W  =   WWWIKK)  »   6. 

281852 

IF(W)  2056  ,2055, 

2056 

2055 

TRRIKK,J,K)  =  0. 
TRKKK, J, K)  =  0. 
GO  TO  2500 

2056 

WTl  =  W«T 
CIS  =   COSF(  WTl) 
SIS  =   SINF(  WTl) 
WDT  =  W  »  TA2 
DCl  =   :OSF{  WDT) 
Dol  =   SINF{  WDT) 

COl  =       1.  /(W«W  « 

TA2) 

WTM  =  W»TAM 

CMS  =   COSF(WTM) 

SMS  =   SINFiWTM)  ■ 

WTL  =  WTM 

CNS  =  COSF  (  WTL  ) 

SNS  =  SINF  (  WTL  ) 

CS  =  1. 

SS  =  0. 

Rl  =  0. 

SI  =  0. 

SH  =  0. 

RH  =  COl  «{  CMS»TLR  +       TCR)  ♦  SMS  *  TLI/W 

IF<  IT)   2060,  2070,  2060 
2060  RH  =  RH  +   CNS  «  TFR  *  COl  +  SNS  »  TFI  /  W 

SH  =  ICMS»TLI  -  CNS«TFI)/W  +  COl» t-SMS»TLR+SNS*TFR ) 
2070  CONTINUE 

M3  =  M  -  2 
:      THIS  LOOP  SWEEPS  THROUGH  CHANGE  IN  SLOPE  CURVE 

DO  2400  KF  =  I,  M3 

CSS  =  CS 

CS  =   CS  •  DCl  -  SS  »  DSl 

SS  =   SS  •  DCl  +  CSS»  DSl 

Rl  =  Rl   +   CS  «  ACJKF,J,K) 

IF{  IT)  2110,  2120,  2110 
2110  CONTINUE 

Rl  =  Rl  +  CS  *AC(KF,K,J) 

SI  =  SI  -  SS»JAC(KF,J,K)  -      AC{KF,K,J)    ) 
2120  CONTINUE 
2400  CONTINUE 

RPS  =  (  RH  +  Rl  »  COl)  /  STRl 

SIP  =  (  SH  +  SI  *  COl)  /  STRl 
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TRR(KK,J,K)  = 

=C1S  » 

RPS  + 

SIP  ♦ 

SIS 

TRHKKtJtK)  = 

=  CIS 

*  SIP  ■ 

-  RPS 

•  SIS 

2500 

CONTINUE 

3000 

CONTINUE 

3050 

CONTINUE 
M  =  NCTl 
NFREQ  =  NWIN 
MZ  =  NTRY 
RETURN 
END 

SPECIAL  NXTPR  FOR  COMBINING  RUNS,   KYLSTRA 

SUBROUTINE  NXTPR 

COMMON  FOR  VARIABLES  SAME  IN  DIFFERENT  MZ  CASES 

COMMON  X,H,G,NCT,ID,ID2,FNT,NRD,NDS,M,MZ,MM,MZ2,MPL 

COMMON  NFILE.WMAX.NDEX, NFREQ, GT 

COMMON  VCTYI.VCTYKtVCTAC 

COMMON  JVAtJVB 

DIMENSION  Xt600) ,H{60I),G(600) 

DIMErjSION    VCTYI(AO),VCTYKUO)  ,VCTAC(40) 

DIMENSION  GT( 18) 

DIMENSION  JVA(40),JVB(40) 

COMMON  FOR  VARIABLES  THAT  DIFFER  FROM  MZ  TO  MZ 

COMMON  Y,AC 

COMMON  TRR.TRI 

DIMENSION  Y(750,4),  AC{300,4,4),  TRR(300,4,4) 

DIMENSION  TRI(300,4,4) 

EQUIVALENCE  {  X,  W) 

DIMENSION   08(300),  R{300),  W(300).  THU300),  TH2(300) 

DIMENSION     T2{300),  TR2(300),  112(300) 

CT  =  57.29578 

CL  =  4.34294482 

NTRY  =  MZ 

NWIN  =  NFREQ 

NTIN  =    2 

NTOU    =  2 

NTPU    =  3 

REWIND  NTIN 

IF  (XABSF(  NWIN)  -  5)  505,  505,  502 

502  DO  503  I  =  1,M 

503  XII)  =  XI n    »  16. 
505  CONTINUE 

DO  50  J  =  1,  NDEX 

K  =  JVA(J) 

L  =  JVB(J) 

DO  521  =  1,  M 

A  =  TRR(  I,  K,L) 

B  =  TRK  I,K,L) 

IF  I  K-L)  4,  5,  4 
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4  CONTINUE 

TH  =  CT»ATANF  I  B/AJ 

THl(I)  =  TH 

TH2( I)  =  180.  -  TH 

AA  =  SQRTFl  A»A  +  B»B) 

Rll)  =  AA 

GO  TO  6 

5  CONTINUE 
AA  =  A 

6  CONTINUE 

XI     =  CL«  LOGF(AA) 
IF{I-1)  2,  2,  3 

2  XS  =  XI 

3  CONTINUE 

DB{  I)  =  XI  -  XS 
52  CONTINUE 

IF(  K-L)  51,  55,  51 
51  CONTINUE 

WRITEOUTPUTTAPE  6,  1001 
1001  FORMAT  {  IHl  lOX  33HVAL 
13HN0RMALIZED  TO      E 
lOX  8H  RUN  NO.  15,  12 
15,  // 
6X  4HFREQ  UX  4HREAL 
lOX  5HTHETA  lOX  5HTHET 
WRITEOUTPUTTAPE  6,  1002 
,  R(  I)  ,  DB(  I  ),  THK  I), 
FORMAT  {  7E15.6  ) 
GO  TO  50 
CONTINUE 

WRITEOUTPUTTAPE  6,  1003 
FORMAT  (  IHl  lOX  31HVAL 
Inh  NORMALIZED  TO    E2 
lOX  8H  RUN  NO.  15,  12 
15,  // 
lOX  4HFREQ  12X  5HRATI0 
5HRATI0  13X  2HDB   ) 
WRITEOUTPUTTAPE  6,  1004 
FORMAT  (  4X  3E17.8,  lOX 
CONTINUE 

IF  (NTRY)  100,  150,  126 
=  -,  WRITE,  =  +,  READ 
CONTINUE 
WRITETAPE  NTOU, 
WRITETAPE  NTOU, 
1  I  =  1,  MJ 
GO  TO  500 
CONTINUE 

READ  TAPE  NTIN,  ID2,  M 
READ  TAPE  NTIN,  {TR2( I) 
NTRY  =  NTRY  -  1 


1002 

55 

1003 


1004 
50 


100 


126 


,  H{ J),  ID,  K,  L 

UES  FOR  CROSS  POWER  SPECTRUM, 

20.8,   // 

H  ,  CHANNELS   15,  7H  ,  AND 

IIX  4HIMAG  lOX  3HAMP  9X  5H  DB 
A    ) 
,  (Wd),  TRR(I,K,L),  TRI(I,K,L) 

TH2( I  ),  I  =  1,  M) 


,  H{ J),  ID,  K,  L 

UES  FOR  AUTO  POWER  SPECTRUM 

0.8,  // 

H  ,  CHANNELS   15,  7H  ,  AND 

13X  2HDB  24X  4HFREQ  12X 

,  (W[  n  ,TRR{  I,;<,L)  ,DB(I)  ,I  =  1,M) 
3E17.8  ) 


ID,  M 

{    TRR{I,1,2),  TRK  1,1,2), TRRd, 1,1), 


TI2(  I),T2( I),  I  =1,m; 
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IFl  NTRY)  150,  150,  127 

127  NTRY  =  NTRY  »  M  «  4 
DO  128  I  =  1,  NTRY 

128  READ  INPUTTAPE  NTPU,  1115,  ITT 
1115  FORMAT  (ID 

150  CONTINUE 

DO  531  =  1,  M 
F3  =  TRR(  I,  1,1  ) 
FR3  =  TRR{ 1,1,2)   /  F3 
FI3  =  TRK  I,  1,2)   /  F3 
F2  =  T2{  I) 
FR2  =  TR2( I)  /F2 
FI2  =  TI2(  I)   /  F2 
DEM  =  FR2»FR2  +  FI2«FI2 
RR  =  (  FR3»FR2  +  FI3»FI2)  /  DEM 
SS  =  (  FI3»FR2  -  FR3»FI2  )/  DEM 
R(  I)  =  RR 
G{ n  =  SS 

AM   =  RR»RR  +  SS»SS 
DB{  I)  =  CL  *  LOGF  t  AM) 
Hd)  =  SQRTF  {  AM) 
THE  =  CT  «  ATANF  (  SS/RR) 
THK  I)  =  THE 
TH2{ I)  =  180.  -  THE 
53  CONTINUE 

WRITEOUTPUTTAPE  6,  1006,  ID  ,  ID2 

1006  FORMAT  (  IHl  20X  30HTRANSFER  FUNCTION  BETWEEN  TWO 
1  lOhDETECTORS   / 

1  lOX  ^2HRATI0  OF  CROSS  POWER  SPECTRUM  FOR  RUN  NO. 
1  15,  22H   TO  THAT  OF  RUN  NO.    15,  //  lOX  4HFREQ  13X 
1  4HREAL  13X  4HIMAG  13X  3HAMP  15X  2HDB  13X  5HTHETA  12X 
1  5HTHETA    ) 

WRITEOUTPUTTAPE  6,  1007,  (W(I),  R(I),  G(I),  H(I), 
1  DB( I),  THK I),  TH2(I),  I  =  1,  M) 

1007  FORMAT  (  IH   7E17.8   ) 
ITT  =  0 

ITH  =  0 

ID3  =  ID«10  +  ID2 

PUNCH  TAPE  B3  FOR  LEAST  SQUARES  INPUT 

DOllO  I  =  1,  M 

WRITEOUiPUTTAPE  NTPU,  1010,  I  TT,  ITH,  W  (I )  ,   R{I),I,Il)3 
110  CONTINUE 

ID3  =  ID3  »  10 

D0115  I  =  :,  M 

WRITEOUTPUTTAPE  NTPU, 10 10 , I TT , I TH, W ( I) ,   S(I),I,ID3 
115  CONTINUE 
1010  FORMAT  (  211,  2E20.8,  2110  ) 

103  =  ID3  «  10 

DO  120  I  =  1,  M 

WRITEOUTPUTTAPE  NTPU , 1010, I TT , I TH, W (  I )  ,  DB{I),I,ID3 
120  CONTINUE 
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ID3  =  ID3  ^    10 

DO  125  I  =  1,  M 

WRITEOUTPUTTAPE  NTPU, 1010, I TT , ITH, W ( I ) , THl ( I ) , I , ID3 
125  CONTINUE 

IF  (  NWIN  )  146,  146,  135 

=  0,  DONOT  REWIND,  =  1,  DO 
135  REWIND  NTPU 

PRINT  1127 
1127  FORMAT  (  31H0MAY  DISMOUNT  TAPE  B3  AND  SAVE   ) 
146  CONTINUE 
500  CONTINUE 

RETURN 

Ei\D 


APPENDIX  E 
COaNT  RATE  CIRCUIT  DESIGN 

The  count  rate  circuit  is  basically  a  simple  diode 
pump  and  storage  capacitor  network,  as  described  in  stan- 
dard textbooks  on  electronic  circuits  {H),     However, 
the  specifications  of  a  fast  response  circuit  and  little 
ripple  on  the  output  signal  cannot  be  satisfied  by  the 
conventional  count  rate  meter  design. 

Conventional  meters  use  large  time  constant  cir- 
cuits to  obtain  an  average  count  rate,  whereas  in  this 
study,  the  ability  of  the  circuit  to  accurately  follow 
the  random  fluctuations  of  the  neutron  density  is  im- 
portant. The  frequency  response  of  the  count  rate  cir- 
cuit is  an  indication  of  how  closely  the  output  signal 
of  the  count  rate  circuit  approximates  the  instantaneous 
variation  of  the  neutron  density.  The  circuit  shown  in 
Figure  32  satisfies  the  specifications,  if  the  design 
criteria  listed  below  are  used  in  selecting  the  values  of 
the  resistors  and  capacitors  in  the  circuit.  References 
(35.  36)  also  provide  information  applicable  to  count 
rate  circuit  design. 
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TABLE  4 
Design  Criteria 


1.  pulse  width  =  1    ^sec 

2.  R^j^  =  lO^ohm 

3.  Rl  =  lO^ohm 

4.  R.^C^  =  10**^360 

5.  (C,+   C2)/(CiC2R,^f^^):>10 

6.  C2(R2+    R3)  =  4/f».^ 


7.     C,(R2+R3)=   1/f* 


max 


8.  C3=    C2/2 

9.  R3/R2  =   2/3 


^°-  ^3  used  =  (^2^  ^3calJ/(^L-^3cal) 
*f  =  counts /sec 

Design  Criteria  1,  2,  and  3  were  established  by  the 
expected  input  and  output  impedance  of  the  equipment  con- 
nected to  the  count  rate  circuit  and  by  the  characteristics 
of  the  input  pulse.  In  practice,  R^^  may  have  any  value 
up  to  10  ohms.  The  circuit  will  continue  to  operate  for 
higher  values  of  R^^i,  but  the  output  voltage  will  be  less 
for  the  same  count  rate.  The  value  of  the  load  resistor, 
Rl#  is  also  flexible.  However,  to  prevent  loading  of  the 
count  rate  circuit,  Rj^  should  not  be  much  less  than  1 
megohm,  particularly  for  the  low  count  rate  ranges.  It 
may  have  any  value  above  1  megohm. 

The  operation  of  the  circuit  is  as  follows.  A  nega- 
tive voltage  pulse  enters  the  circuit  and  starts  charging  C, 
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As  C,  is  charging,  the  arrangement  of  the  diodes  is  such 
that  T)^   shorts  C,  to  ground.  The  equivalent  circuit  is 
a  series  circuit  of  R^^^  and  C  .  Thus,  if  Design  Criterion 
4  is  satisfied,  C,  is  completely  charged  by  the  one  micro- 
second long  pulse  supplied  by  the  linear  amplifier.  When 
the  pulse  stops,  C]^  discharges  through  D2,  transferring 
part  of  its  charge  to  Cg.  If  Design  Criterion  5  is  satis- 
fied, this  transfer  is  easily  made  before  the  next  pulse 
occurs.  Although  Gg  continuously  discharges  through  ^2'^"^^ 
Rg,  Design  Criterion  6  prevents  complete  discharge  of  Cg 
before  the  next  pulse  occurs. 

Design  Criterion  7  is  used  to  insure  that  the  output 
of  the  circuit  is  approximately  proportional  to  the  input 
pulse  rate.  The  output  voltage  is  generally  less  than 
25  per  cent  of  the  peak  voltage  of  the  incoming  pulses. 

If  C3  is  not  in  the  circuit,  the  output  signal  fluc- 
tuates considerably  between  pulses.  This  fluctuation  is 
caused  by  the  fast  response  of  Cg  and  the  bypass  resistors 
R*  and  R-.  The  presence  of  C«  adds  a  low  pass  filter  that 
smooths  the  ripple  without  reducing  the  frequency  response 
of  the  circuit.  Choosing  the  values  of  Ti  ,   R  ,  and  Cg 
according  to  Design  Criteria  8  and  9  will  reduce  the  ripple 
on  the  signal  from  the  basic  circuit  by  more  than  95  per 
cent  and  increase  the  frequency  response  of  the  circuit 
one  and  one  half  to  two  times.  Design  Criterion  10  adjusts 
the  value  of  R3  to  condensate  for  the  loading  effect  of  Rj^. 

Switch  1  has  two  rotary  wafers,  allowing  simultaneous 
changing  of  R2  and  R3  to  obtain  a  different  count  rate 
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range.  If  desired,  C^  may  be  included  in  the  circuit  by- 
Switch  2  to  remove  the  dc  level  from  the  output  signal. 
The  values  of  the  components  used  in  the  count  rate  cir- 
cuit are  tabulated  in  Table  5. 


TABLE  5 

Count  Rat' 

e  Circuit 

Components 

Range 

Counts  per  sec 

24  K 

R3,  ohm 
16  K 

-3db  point 

1 

100  K  to  250  K 

5400  cps 

2 

40  K  to  100  K 

60  K 

42  K 

2000  cps 

3 

20  K  to  50  K 

120  K 

87  K 

1000  cps 

4 

10  K  to  25  K 

240  K 

195  K 

600  cps 

5 

4  K  to  10  K 

600  K 

667  K 

350  cps 

Ci  =  100  p^pif  Cg  =  1000  /i/if 

C3  =  500  pLpLf  C^  =    2  pipLf 

D^  and  D2  are  high  speed  silicon  computer  diodes. 

Figure  33  shows  the  attenuation  of  the  circuit  when  on 
Ranges  1,  3,  and  5. 
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